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Abstract 


Two approaches are developed to analyze the dynamic behavior of flexible multibody 
systems. In the first approach each body is modeled with a modal methodology in 
a local non-inertial frame of reference, whereas in the second approach, each body is 
modeled with a finite element methodology in the inertial frame. In both cases, the 
interaction among the various elastic bodies is represented by constraint equations. The 
two approaches have been compared for accuracy and efficiency: the first approach is 
preferable when the nonlinearities are not too strong but it becomes cumbersome and 
expensive to use when many modes must be used. The second approach is more general 
and easier to implement but could result in high computation cost for large system. 
The constraints should be enforced in a time derivative fashion for better accuracy and 
stability. 
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Chapter 1 
Introduction. 

A comprehensive analysis methodology for dynamic systems involving several elastic 
bodies must include a scheme to efficiently deal with the interaction between the various 
elastic components. For instance, in a conventional helicopter, the elastic rotor or rotors 
interact with the elastic fuselage, whereas in a tilt rotor configuration, the elastic rotors 
interact with the flexible wings and fuselage. The impact of this interaction is important to 
an accurate prediction of rotor loads, and essential when attempting to predict instabilities 
such as ground or air resonances since rigidly mounted rotors do not exhibit such 
instabilities. It is convenient to think of the helicopter as a multibody elastic system, 
i.e. a collection of elastic bodies mutually interacting at “hinges”. 

A fundamental difficulty in the analysis of a multibody system is the evaluation of its 
total kinetic energy, as it involves the calculation of the inertial velocity of each material 
point of the system. If the position of all material points is measured in a given inertial 
system, this task is trivial, however, it is often convenient to use a local coordinate system 
to represent the initial geometry and deformation of each elastic body. The velocity of 
a material point relative to this local frame is easily to obtain within the frame work of 
a finite element discretization or modal representation, however, the inertial velocity of 
this material point also involves the motion of the local frame with respect to an inertial 
frame of reference. This additional motion can be taken into account through various 
schemes, for instance hierarchical representations, or multibody schemes. 

A hierarchical representation involves a hierarchy of reference flames starting with 
an inertial frame. The motion of each frame is described with respect to the frame that is 
immediately superior in the hierarchy. For a helicopter, a typical hierarchy could be as 
follows: inertial frame — to — airframe system — to — blade system — to — deformed 
blade system. Each level of the hierarchy involves a rotation matrix which gives the 
instantaneous position of a frame with respect to that immediately superior. Each rotation 
matrix is quadratic in terms of the Euler Parameters (other finite rotations parameters 
could be used but the rotation matrix will remain at least quadratic). Since our typical 
hierarchy involves four levels, the position and inertial velocity vectors of a material 
particle will involve nonlinear terms up to the 9th order, resulting in a kinetic energy 
expression with nonlinear terms up to the 18th order. Of course, some simplification could 
be introduced: for instance in level flight, the inertial — to — airframe transformation 
becomes a constant, or, if the elastic deformations of the blade are linearized, the blade 
— to — deformed blade transformation becomes linear. However, a general analysis 
methodology should be able to deal with large rotations at all levels. 

This simple example points out the two major difficulties associated with hierarchical 
models: first, these models are difficult to handle and require advanced data base concepts 
for practical implementation, and second, very high order nonlinear terms appear in 
the analysis resulting in a very large number of coefficients. In a modal analysis, the 
number of coefficients is N n , where N is the number of modes, and n the power of the 
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nonlinearities. For a 12 mode model involving 18th order nonlinearities 2.66 10 19 terms 
will be generated, requiring 2.03 10 14 Mbytes of storage on a computer. Of course the 
number of operations involved in manipulating the model grown accordingly. It is clear 
that such model are beyond the reach of even the most powerful computers, and would 
require exorbitant amounts of computation. 

Two alternative approaches will be pursued in this work that avoid hierarchical 
representations. In both approaches, the dynamic system is modeled as a collection of 
flexible bodies ( airframe, wings, blades, etc...) that are connected together at a number 
of points where kinematic constraints are enforced. Typical kinematic constraints are 
spherical, universal and convolute joints, or rigid links. In the first approach each elastic 
body is described in a local coordinate system which motion is directly related to an 
inertial frame through three rigid body translations and three rigid body rotations. Hence, 
the inertial position vector of any material particle in that elastic body involves a single 
rotation matrix only, allowing an easy evaluation of all inertia terms. 

Since a local coordinate system is used, the elastic deformations of the body can 
be represented in a modal fashion, more specifically a finite element based modal 
analysis technique will be used which yields a Lagrangian expression involving quartic 
nonlinearities only. 

In the second approach, all elastic bodies are described directly in a single inertial 
system. This is by far the simplest formulation, however, it rules out the use of a modal 
representation, and requires a parametrization of the finite rotation variables that allows 
arbitrarily large rotations ( in this work the Milenkovic parameters are used.) 

Chapter 2 presents a review of the beam model which will be used throughout this 
work. The next two chapters deals with the modeling of a single elastic body: Chapter 3 
presents the modal reduction scheme for the finite element model, and chapter 4 compares 
the predictions of modal models with that of full finite element models. The kinematic 
constraints to be applied between elastic bodies is the focus of chapter 5. Chapter 6 briefly 
describes the full finite element modeling. Finally conclusions and recommendations are 
presented in chapter 7. 
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Chapter2 


Finite Element Modeling of Rotor Blades 
Section 2.1: Introduction 

The kinematics involved in the nonlinear static and dynamic analysis of naturally 
curved and twisted blades are complex since both the deformed and undeformed config- 
urations of a blade are three dimensional. Moreover, laminated composite materials are 
increasingly used for the construction of such structures, causing several non-classical 
effects of beam theory to become more pronounced [2-1,2]. 

In many applications, large displacements and rotations of the blade will occur; 
however, the strain level remains low. Fatigue life is indeed a major concern; hence, the 
operating strain level must remain well within the linear-elastic range of the material. As 
a result, most analyses [2-3,8] are based on a small strain assumption that considerably 
simplifies the formulation and resulting equations. 

The small strain assumption has important implications. First, the Green-Lagrange 
strain components often used in the derivation of nonlinear kinematics [2—4,5] can 
be equated to the engineering strain components, and hence the usual stress-strain 
relationships of the material can be used. Second, the changes in surface area of 
a differential volume element due to deformation arc neglible. Finally, the strain- 
displacement equations can be considerably simplified, since all second order terms (i.e., 
strain square terms) can be neglected. 

In this chapter, consistent strain-displacement expressions are derived which provide 
the basis of the finite element approximation of the non-linear behavior of naturally 
curved and twisted blades undergoing arbitrarily large deflections and rotations. 

Section 2.2: Geometry and Kinematics of Blade Elements 

Consider a naturally curved and twisted beam depicted in Fig. 2.1.1 The triad Ti , i 2 , 13 
is fixed in space and the triad ej, e*2, e 3 is attached to a reference line along the axis of 
the beam. t\ is chosen tangent to the reference line and e 2, e 3 define the plane of the 
cross-section. The curvilinear coordinates along this triad are xi,x2, 23 respectively. 

The position vector of a particle of the beam in the undeformed configuration is: 

f = r (xj, x 2 , £3) (2.2.1) 

After deformation the same particle has a position vector: 

R = R(x i,x 2 ,x 3 ) (2.2.2) 
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Reference Line 



Figure 2.1.1. Geometry of the Beam Before and After Deformation 
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The corresponding vectors at the reference line are: 

1 = * (2.2.3) 

Ro = i?o (^i) 0, 0) 

and the displacement vector of the reference line is given by: 

u = Ro - fo (2.2.4) 

The base vector [2-9,10] in the undeformed and deformed positions respectively are 
defined as: 


g t = r*i and G, = if, (2.2.5) 

where the notation (.) ■ means derivative with respect to x,. At the reference line the 
base vectors are: 


e*i = r 0 7, and E t = R 0> i 


( 2 . 2 . 6 ) 


e, forms a triad since the derivatives in (2.2.6) are taken with respect to the natural 
coordinates of the beam. The triad e, can be viewed as a rotation of the basic reference 
triad i, through a given rotation matrix t T (x i) such that : 


e i 


_ ^ _ 
h 

el 

= t T (xi) 

h 

. e~3 _ 


.*3. 


The derivatives of this triad are readily calculated as : 
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(2.2.7) 


( 2 . 2 . 8 ) 


where the notation ()' means derivative with respect to xi; k 1 is the natural twist (or 
pre-twist), k 2 and fc 3 are the natural curvatures (or pre-bends) of the beam. The position 
vector f of an arbitrary point of the beam can now be written as: 


r = ro + X 2 e 2 + X 3 e 3 (2.2.9) 

hence the base vectors become: 

gl = \fgz i — x 3 &ie2 + X 2 fcie$ 

92 = el (2.2.10) 

93 = el 


where 

y/g = 1 — X 2 kz + X 3 ^2 (2.2.11) 

The metric tensor is obtained as g tJ = g x ■ gj and its determinant is g. 
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The fundamental assumption in beam theory is that the cross-section does not deform 
in its own plane. This means that the base vector E 2 and E 3 which are in the plane of the 
cross-section after deformation simply correspond to a translation and rotation of the base 
vector e *2 and e *3 of the original configuration. Note that arbitrarily large displacements 
and rotations can occur but no deformation of the cross-section is allowed i.e. E 2 and 
E 3 are mutually orthogonal unit vectors. In contrast, E\ is no longer a unit vector nor 
is it orthogonal to E 2 and E 3 , as axial and shearing strains are allowed. Now a new 
orthogonal triad e* is defined as follows: 

e 2 = -^2 

e] = E 3 (2.2.12) 

e * = e 2 x e 3 

The vector E\ can be resolved in this triad as: 

£1 = (1 + en) e* + 2ei 2 e* + 2e 13 e* (2.2.13) 


At this point, eif, en are the unknowns, and they will be identified later as 
strain quantities. Here again the triad e* can be related to the basic reference triad i, 
through an unknown rotation matrix T e (x\) such that: 
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(2.2.14) 


The derivatives of this triad are: 
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(2.2.15) 


where K\ is the twist, K 2 and K$ are the curvatures of the deformed beam. Since the 
cross-section does not deform in its own plane, the position vector R in the deformed 
configuration can be written as: 


R = Rq + x 2 E 2 + x 3 E 3 + 8 (xi) ipi (x 2 , x 3 ) ej (2.2.16) 

The first three terms represent large translations and rotations of the cross-section and 
can be geometrically interpreted as plane sections remaining plane, but not necessarily 
normal to deformed axis of the beam(i.e. a Timoshenko Beam Theory). The last term 
represents a small displacement in the direction of e], that is out of plane warping of 
the cross-section chosen as the torsion related warping displacement tpi. This warping 
displacement is selected as the Saint- Venant torsional warping functions[2-l]. 6 (ij) is 
an unknown function characterizing the magnitude of the torsional warping. Combining 
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(2.2.5), (2.2.6), (2.2.13) and (2.2.15), the base vectors of the deformed configuration 
become: 

G\ = [(1 + en) — X2K3 + X3K2 + *V,] ^1 + [ 2 c 12 — X 3 A'i] 

+ [2ei3 + X2K\] &S (2 2 17) 

G2 = ^ 1 , 2^1 + <?2 
Gz = 8y>\,zc\ + 

In (2.2.17), all higher terms containing warping quantities have been neglected. 


Section 2.3: Strain Analysis 


The Green-Lagrange strains f X} in the curvilinear coordinate system [2-9,10] are 
given as /, j = | (G X] - g x j ) where G x} = G x • Gj is the metric tensor in the deformed 
configuration. It is straightforward to verify that f 22 = /33 = /23 = 0 as a direct 
implication of the indeformability of the cross-section in its own plane. The other strain 
compoments are the two transverse shearing strains /12 and /13 , and the axial strain fn- 
To relate these strains to the strains in the local rectangular coordinate system defined 
on the beam axis, the following transformation is needed. Define a local rectangular 
cartesian coordinate system y x along e„ then the relation of this rectangular system with 
the material coordinate system xj is governed by: 


dy x 

dxj 


' y/g 0 0 

— X3&1 1 0 

X2&1 0 1 


(2.3.1) 


Now, the strains e XJ defined in the local rectangular coordinate system e, are obtained as: 


-tj 


dik dxi 
dyidyj 


fki 


Then, the non-vanishing strain components become: 


\/ff e 12 = /12 

y/<j e 13 = /l3 

\[gz 11 = /ll + 2X3 ^/! 2 - 2x 2 fci/i3 


(2.3.2) 


(2.3.3) 


The initial curvatures of the beam ^ and k$ are now assumed to be small, i.e. from 

(2.2.11): 

y/g « 1 (2.3.4) 

In the case of helicopter blade, 


X 2 ^3 


X3k 2 


chord of the blade 

radius of the edgewise curvature 
thickness of the blade 

radius of the flapwise curvature 
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Note that y/g = 1 for a straight blade. Hence, this assumption is realistic for most of 
the practical applications. 

The strain components now become: 


2 ei 2 = 2 ei 2 - X 3 /C 1 + 5^1,2 

(2.3.5) 

2ei3 = 2ei2 + X2«i + ^ 1,3 

(2.3.6) 


e il — cn + -efj — X 2 (1 + en) «3 + X3 (1 + en) «2 + ^Vl + ^1 (^3^1,2 — ^2^1,3) 
+ 2 (2C12 — ^3«l ) 2 + 2 (2ei3 + x 2 «l ) 2 + 2 ( x 2&3 — X 3 /C 2) 2 

(2.3.7) 

where K t = K t - ki 

To complete the formulation, the coefficients en, efi, in (2.2.13) must now be 
related to the displacements and rotations. Differentiating (2.2.4) with respect to xj and 
using ( 2 . 2 . 6 ), we obtain: 


E\ — e\ + u 


(2.3.8) 


or 

El = (u'j + <n) il + (u ' 2 + < 21 ) *2 + (uj + <31 ) Q (2.3.9) 

where u, are the components of the displacement vector in the basic reference triad i t , 
and tij the components of the rotation matrix t. On the other hand, combining (2.2.13) 
and (2.2.14) yields another relation for E\ that can be identified with (2.3.9) to obtain 


1 + en 
2ei2 
2ei3 



uj + <n 
+ <21 
u 3 + <31 


(2.3.10) 


This completes the strain analysis. It is important to note that this development is valid 
for arbitrarily large displacements, rotations and strains. The unknowns of the problem 
are the three displacements u,, the rotation parameters implicitly defined in the rotation 
matrix T e , and torsional warping amplitude. 

In the derivation of strain expressions in (2.3.5— 7), no assumptions were made about 
the magnitudes of the displacements, rotations or strains; hence, these expressions are 
valid for beams with small initial curvatures undergoing arbitrarily large displacements, 
rotations and strains. On the other hand, later in the derivation of the total potential 
energy expression in (2.4.17), strain components will be assumed small enough to render 
negligible changes in area due to deformation, and to equate the Green-Lagrange strains 
to engineering strains. This requires both axial and shearing strains to be much smaller 
than unity, i.e. e« 1 , and 7 « 1 . However, nothing is assumed about the relative 
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magnitude of e versus 7 . For consistency, the same assumptions must now be introduced 
to the strain-displacement equations (2.3.5-7) to obtain: 

£11 =en - x 2 k 3 + z 3 « 2 + S'tpi + Ski (x 3 ip lt2 - x 2 </?i )3 ) 


+ 2 (^ e i 2 - x 3 K i ) 2 + 2 ( 2 ei 3 + X 2 «i ) 2 

(2.3.11) 

712 = 2ei2 - x 3 k\ + fy?i >2 

(2.3.12) 

713 = 2 ei 3 + x 2 k\ + fy?i j3 

(2.3.13) 


The last term appearing in (2.3.7) is negligible since it represents the square of the axial 
components due to bending. 

If we now introduce the additional assumption that axial and shearing strains are 
of the same order of magnitude, then j 2 « e, and the two last terms in ( 2 . 3 . 11 ) can 
be neglected, since they are squares of the shearing strain components in ( 2 . 3 . 12 ) and 
(2.3.13), respectively; this yields: 

£11 = en - *2*3 + *3*2 + + Ski (x 3 <pi t 2 - X2W\,z) 

712 = 2ei2 — x 3 «i + 6^1,2 (2.3.14) 

713 = 2 ei 3 + x 2 k\ + S<fi i3 

These expressions are often successfully used as the basis for beam models involving 
large displacements and rotations, but small strains [2-4 to 8 , 2-11,12]. However, it is 
interesting to note that one additional assumption was required (y 2 « e), that might not 
be adequate when dealing with highly anisotropic composite materials [ 2 - 12 ], 


Section 2.4: Blade Strain Energy Expression (Hellinger-Reissner Formulation) 


The strain energy expression for a thin walled beam is: 


L 

U = ^ J J sLAzdsdxi (2.4.1) 

0 r 

where £ = (e, 7 ); A is the stiffness matrix; L the length of the blade; T the contour of the 
thin-walled section, described by a curvilinear variable s (see Figure 2 . 4 . 1 ). Consistent 
with the assumption of a cross section that does not deform in its own plane, the only 
non-vanishing strain components are the axial strain e and the shearing strain 7 . Clearly, 

£ = £11 and 7 = x 2 * 7 x 2 + x^ 7 J3 (2.4.2) 

where () + denotes a derivative with respective to s. A set of strains e = £ is now 
introduced into (2.4.1) to yied: 


U 


0 r 


,r at 


A 1 e — n_ (e — e) 


ds dx 1 


(2.4.3) 
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where the condition e = £ was enforced by means of a set of Lagrange Multipliers n. 
Variation over e yields n = Ae which yields the physical meaning of n as the internal 
stress flows n = ( n,q ) , where n is the axial stress flow and q the shear stress flow. 
Finally the strains e arc eliminated from (2.4.3) using e = A* 1 n to find: 


-// 


„T _ 1 T A -\ „ 

n e n A n 

“ 2 — 


ds dx i 


Introducing (2.3.14) into (2.4.2) yields 


(2.4.4) 


£ = e n - x 2 k 3 + x 3 k 2 + tpi# — kiripfS ^ 

7 = x 2"2ei2 + £j"2ei3 + r/ci + 

where r = x 2 x^ — x 3 x^ is the distance from the origin to the tangent to the contour r 
(see Figure 2.4.1); tpf = x^v?i, 2 + xf<p i >3 and x 3 y?i, 2 - x 2 ¥>i ,3 - 



Figure 2.4.1 Geometry of the Thin- walled Cross-section 


Note that equation (2.4.5) imply the small strain assumption and can be written as: 


£ 


1 0 

0 

-kir<pf 0 i3 — x 2 (fix 

7 


.0 4 

4 

4 r 0 0 0 


(2.4.6) 


where 

EL- — (ci, e 2 , e 3 , e 4 , ki, /c 2 , « 3 , * 4 ) 

and (2.4.7) 

e l = cii, e 2 = 2ei 2 , e 3 = 2ei3, e 4 = S, /c 4 = <$' 
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Introducing these equations into (2.4.4) and integrating over the cross-section yields: 

L 


U = J [eiFi + C.2F2 + ezFz + e 4F4 + k\M\ + K2M2 + /C3M3 + K^Mildxi 


Lt 

-// 

0 r 


(2.4.8) 


1 I „T A - 1 „ 

— i n A n 

2 


ds dx 


1 


where 


F\ = J nds, F2 = J x^qds, £3 = J x^qds , 
r r r 

F4 = J pf (q — kirn) ds , Mi = J rqds , 

r 

— J X2nds, M4 = j (finds 


(2.4.9) 


r r 

M2 = / x^nds, M3 

r r r 

Fi is the axial force, Mi is the torque, F2 and F3 are the shear forces, M2 and M3 are 
the bending moments, and finally F4 and M4 are the force and moment, respectively, 
associated with the torsional warping induced stresses. 

With Reissner’s Principle independent assumptions can be made on displacements and 
stresses. By analogy to the strain field the stress field is assumed in the following form 
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^nn 

Anq 
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. ^nq 

■^99 . 



L ^"9 ^99 J 
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B X 
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0 — kir<fi 0 X3 — X2 ip 1 

• + -+ r 0 0 0 


<Pi 


X 


(2.4.10) 


where X T = {Xi , X 2 , Xz, X4, Yi , Y2, F3, F4) is a vector of unknown stress parameters, 
P 2 and fz the transverse shearing related Saint- Venant warping functions[2-l]. Intro- 
ducing (2.4.10) into (2.4.9) and integrating over the cross section yields 


£ = AX 


(2.4.11) 


where EL = (£1 , £2 , £3 , £4 , Mi , M2 , M3 , M4 ) and A is a matrix of cross-sectional 
coefficients. Finally the stress assumptions (2.4.10) are introduced into (2.4.8) to find: 


it 

U = J [(e,Fi + e2£2 ■+• C3F3 + 64^4 -j- kiMi + /C2M2 + K3M3 4 - K4M4) 


(2.4.12) 


- \ ELH E] dx\ 

where the compliance matrix is H given by 

H = A~ T DA- 1 


(2.4.13) 
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where 


D = 




B ds 


Finally, (2.4.12) can be written as: 


L 

v = J (£e.-\eLh A dx i 

0 ' 


(2.4.14) 


(2.4.15) 


where e£ = (ei,e2,e3,e4,«],K2»«3j«4) arc the sectional strains, F£ = 

(Fi , F 2 , F 3 , F 4 , M 2 , M 3 , M 4 , F 5 ) are internal forces, and the strain-displacement re- 

lationships are: 


'1 + ef 


<11 +«i 
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(2.4.16) 
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This nonlinear strain energy expression depends on the displacements and forces: 


U = U (n,E). 

It can be expanded using a linearization procedure about a reference configuration, to 
yield: 


where 


and 



Uuu 


<¥_U 
dudu ’ 


U uf = 


d 2 U 

dudF 


U ff = 


d 2 U 

dFdE 


(2.4.17) 

(2.4.18) 

(2.4.19) 
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Section 2.5: Blade Kinetic Energy Expression (Hellinger-Reissner Formulation) 



Let Ii be an inertial reference frame. Consider now an unstrained structure in space 

with a triad i t attached at a material point O. To locate the structure, it is convenient to 

separate its displacement field into rigid body displacements and elastic displacements. 

The rigid body motions define the position of a fictitious, rigid structure, and the elastic 

motions are superimposed to yield the true position of the structure. The rigid body 

displacement field involves three translations and three rigid body rotations. The rigid 

body translations are chosen as the translations of a material point O with an unknown 
— # 

position vector P ( t ). The rigid body rotations consist of two parts: first, an unknown 
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rigid body rotation characterized by a rotation matrix T r (<), then a known rigid body 
rotation with constant angular velocity charaterized by a rotation matrix T r 0 (t) such that: 


Jl ' 

32 

= 2f (0 

'/f 
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’{l " 
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ii 
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.33. 
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which yields: 
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= T? T {t)Tj{t) 
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. *3 . 
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(2.5.2) 


The elastic displacement field involves elastic displacements u{x\,t ) of the reference 
line defined in (2.2.4) and elastic rotations as defined in (2.2.14). 

The sum of rigid and elastic displacement fields brings the structure to its actual 
position. The fictitious, rigid structure is used as a reference configuration for the 
described of the elastic strain field, in a manner identical to that described in section 2.2. 
This involves a set of material coordinates xi,x 2 ,x 2 , the base vectors of the reference 
line in the undeformed configuration e\, e 2 , e 3 , and the base vectors of the reference 
line in the deformed configuration ej, e^, e^. 

The position vector of a material point is: 

R = P + fo + u + x 2 e 2 * + x 3 e* 3 (2.5.3) 


The instantaneous position of a material point on the blade is given by (2.2.16) where 
torsional warping related out-of-plane displacements were included, however the inertia 
forces associated with this out-of-plane motion are very small and will be neglected here, 
resulting in the simplified expression (2.5.3). It is clear that rigid body motions do not 
generate any strain field, hence the strain displacement equations (2.3.14), and the strain 
energy expression (2.4.12) remain unaffected. However, the kinetic energy expression is 
of course affected by the presence of time dependent, rigid body motions. 

The time derivative, noted (), of (2.5.2) yields: 
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(2.5.4) 


The time derivative of (2.2.14) yields: 


L e 3 J 


= (if, sfr, + ifafr, + if: I-,) 
= (o; r +«•**• +*;!•) 


L e 3 J 


s 3 J 


(2.5.5) 
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The following skew symmetric angular velocity matrices have been defined: u T and 
the unknown and constant rigid angular velocities, respectively, resolved in the i t system; 
u)*, u)j?* and u>* the unknown, constant, and elastic angular velocities, respectively, 
resolved in the e* system. The components of the corresponding angular velocity vectors 
can also be defined: uv = (w r i,w r 2 ,w r3 ), wj? = (a^u^u^); u£ = (u^w^,^), 
= (“>?!,<«$, w°j) , and finally It is clear that: 

u* = Tj ov and = Tj wj (2.5.6) 


The position vector (2.5.3) can be resolved as follows: 


R=[h h h) 


Pi 

P 2 

P 3 


+ *1 *2 h . 


Xo + Ul 

yo + «2 

ZO + U 3 


+ [e] el el] 


0 

x 2 

X 3 


(2.5.7) 


where Pj_ = (Pi,P 2 ,P 3 ) are the components of the vector P in the inertial system. 
The inertial velocity of a material particle is found by combining (2.5.7) with (2.5.4) 
and (2.5.5), to find: 


J2=[*I *2 iz){T? T Tj 


Pi 

P 2 

P 3 


+ (fir + ) 


Xo + «1 


«l" 


yo + U 2 

+ 

U 2 

> 

. 20 + U 3 _ 


_u 3 _ 

J 


+ [ei* er «r ](«;+«;• +«:) 


0 

x 2 

X3 


(2.5.8) 


The following notations are introduced. First, the rigid translational velocities are 
defined as: 



Vrl 

11 

H 

‘Pi' 

El ~ 

Vr2 

P 2 




P 3 


+ X T U r 


(2.5.9) 


where 



" v el " 


'ui ' 


v e2 

= 

«2 


_v e3 _ 


.W3. 


+ X T u>° 


X = 


0 - (20 + U 3 ) (yo + u 2 ) 

(20 + u 3 ) 0 — (xo + «i) 

-(yO + U 2 ) (xo + Ul) 0 


(2.5.10) 


(2.5.11) 


The total translational velocity is now defined in the Ti system as: 


Vt = V r + V e 


(2.5.12) 
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and in the e? as: 


vl = Tjv, 

Finally, the total angular velocity vector in the e* sytem is: 

* * , 0+ , * 
u>t = u T + u r +U) e 


(2.5.13) 

(2.5.14) 


For (2.5.9) to (2.5.12), the inertial velocity of a material point (2.5.8) becomes: 

v t\ ~ x 2 u t3 + x 3 u t2 


R = [eV eT eV] 

The kinetic energy of the system is: 

L 

1 

T 




X 3 U> t *j 


V *Z + *2^*1 


IJJok-i 

o r 


Rdsdx i 


Introducing (2.5.12) and integrating over the cross section yields: 

L 

T = \ j V t * T MV t * dx\ 
o 

where the array of total velocities is defined as: 

Yl _ = ( V *l’ V t2 ' v t3 1 u (] > u t2 1 u t3 ) 

and the mass matrix M is given by: 


M = 


m 

0 

0 

0 

m3 


0 

0 

0 

m 

0 

—m 3 

0 

m 

m 2 

-m 3 

m 2 

mu 

0 

0 

0 

0 

0 

0 


m 3 —m 2 

0 0 

0 0 

0 0 

77733 77723 

— 77723 77722 


where the mass per unit span is: 


m 


-/ 


pds 


the inbalance per unit span is: 


m>i — f px j 


ds i = 2, 3 


(2.5.15) 


(2.5.16) 


(2.5.17) 


(2.5.18) 


(2.5.19) 


(2.5.20) 


(2.5.20a) 
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the moment of inertia per unit span is: 


mij = 


^ p x j i j d s 

r 


i,j = 2,3 


and 


mn = m22 + m33 


(2.5.20b) 


(2.5.20c) 


For future reference, it is convenient to write: 


v; = av l +v; 


where the matrix A is defined as: 


A = 


Tj TjX 7 
0 Tj 


the rigid velocities are: 


T 0T T t p 


T? t 0 

Tj P 

Ur. 


0 T r or 

2GQ 


and the elastic velocities are: 



(2.5.21) 


(2.5.22) 


(2.5.23) 


(2.5.24) 


A set of velocities JJ_ = VS is now introduced into the kinetic energy expression 
(5.2.14) to yield ~ 


L 

T = J \uLMTL- p ' t (H-V*) dx 


(2.5.25) 


where the condition U — V* = 0 was enforced by means of a set Lagrange multipliers 
P* . Variation over yields p*_ = MU_ which yields the physical meaning of p* as the 
momenta components, measured in the system. Finally, the velocities JJ_ are eliminated 
from (2.5.18) using H = M~ l to find: 


T = 



dx i 


(2.5.26) 


this expression of the kinetic energy is a nonlinear function of the six rigid body velocities 
Vt — { v t\,v T 2 , v T z,u T i,u r 2 ,u T {), the elastic displacements u and their time derivatives 
u, and the compontents p* i.e.: 


T = T(K,u,u,p*) 


(2.5.27) 
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The kinetic energy will be used in two way. First, a finite element implementation 
of the problem to yield steady equilibrium position of the structure under applied loads 
and normal vibration modes. Second a modal approximation to obtain modal equations 
of motion. In the finite element implementation which focuses on natural vibration mode 
calculations, constant, rigid angular velocities are the only allowed type of rigid body 
motion. On the other hand, in the modal approximation the 6 rigid body motions are 
unknowns of the problem and describe the rigid body response of the system to the 
applied load. For constant rigid angular velocities we have 


T = T («,«,£) 


(2.5.28) 


this nonlinear expression can be expanded using a linearization procedure. This expansion 
is performed about a steady configuration noted p*° , and the rigid angular velocities 
are to find: 


/{ 


T = l<r + [A« A u A p* ] 
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l T pj 
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— dudiL ’ 
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T uu T u £ T U p 
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Tuu Tuu Tu P 
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_ T U p T^p Tpp _ 
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_ dT 
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— — Pt » 

on 

di A 
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+ H.O.T.^dx 1 


_ d 2 T 
d 2 T 


T pp — 


d£_dp*_ 


All these arrays are evaluated in the reference configuration. 


(2.5.29) 


dT 

7T7 (2-5.30) 


(2.5.31) 


Section 2.6: Normality Condition for Euler Parameters 


In the two previous sections finite rotations were used, and to keep the formulation 
general, the rotation matrix T e only appears in the equations. However, for a practical 
implementation, a specific set of rotation parameters must be selected. In this work 
the Euler Parameters (see appendixA) are used that are related through the normality 
condition. This normality conditon could be enforced using a penalty method, i.e. adding 
the following term to strain energy 


C = 



( 2 . 6 . 1 ) 
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where a is a large penalty coefficient, and 


e 9 — 9o + 9l + 92 + 93 — 1 


( 2 . 6 . 2 ) 


is the normality condition and q, the Euler Parameters. Following the Hellinger-Reissner 
approach used in the previous sections, a variable e = egis introduced into (2.6.1) to 
yield 





- Me - e 9) 


dx i 


(2.6.3) 


where the condition e = eg was enforced by means of a Lagrange Multiplier A. Variation 
over e yields A = ae, hence, the variable e can be eliminated using e = A /a, to yield 


L 



0 


(2.6.4) 


It is convenient to interpret this relationship in the following physical terms : eg is a 
fictitious strain, a is fictitious stiffness, and A a fictitious force, such that A = aeg. 
By selecting a very large stiffness a we drive the strain eg to zero, i.e. we verify the 
normality condition. 

The nonlinear constraint expression (2.6.4) can be expanded using a linearization 
technique to yield 


where 


C = C + [Ai£ AAl] 

+ |[A i£ AAl] 


c A 


Ci tu Atx 

C u \ C\\ AA 


c _££ 

^ u — ^ 

— au 


c, LU = 


d 2 C 


r =?£ 

— dx 

C u \ = 


d 2 C 


■f h.o.t. 


uu dudu ' ~ UA - dud A’ 

All derivatives are evaluated in the reference configuration. 


Cxx = 


d 2 C 
d\d\ 


(2.6.5) 


( 2 . 6 . 6 ) 
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Chapter 3 

Modal Reduction of the Finite Element Model 
Section 3.1: Introduction 

The versatility and efficiency of the finite element method makes it an attractive tool 
for the analysis of helicopter rotor blades. Hodges [3-1] has recently reviewed various 
finite element approaches, giving a comprehensive discussion of their assumptions and 
features. An analysis including moderate rotations was developed by Friedmann and 
Straub [2-3], as by Sivaneri and Chopra [3-2] based on the formulation of Hodges and 
Dowell [3-3]. Giavotto et al. [3—4] and Born [3-5] developed an approach that includes 
finite rotations, as well as cross sectional warping deformations. Finally, a model for 
arbitrarily large displacements and rotations of naturally curved and twisted blades was 
developed in chapter 2. 

These various approaches are very attractive because they allow accurate modeling of 
rotor blades. The complex kinematics resulting from the large displacements and rotations 
can be handled in a rational manner, and the intricate elastic behavior of composite blades 
can be treated realistically by introducing transverse shearing and warping deformations, 
as well as elastic coupling. However, the cost of such analysis can be prohibitive when 
realistic problem must be treated. 

Consider a composite blade with varying properties along the span: 100 to 150 
degrees of freedom (DOFs) are typically necessary to accurately represent its geometry 
and physical properties. This number must appear small, as problems involving 1,000, 
or even 10,000 DOFs are routinely solved with large finite element codes. However, 
in the case of a helicopter blade, the analyst is interested in determining its nonlinear 
static behavior, its dynamic characteristics i.e. its natural frequencies and mode shapes, 
its nonlinear, periodic dynamic response, and the stability characteristics of this periodic 
response. The first two analysis types are relatively straightforward to handle, but the 
latter ones are far more complex. 

Consider the prediction of the nonlinear periodic response of the blade using the finite 
element in the method [3-6]: the total number of DOFs equals the number of DOFs used 
for the spatial model times the number of time stations. If 64 time stations are used, this 
will yield 6,400 to 9,600 DOFs to be solved for in an iterative manner, since the problem 
is nonlinear. For a gimballed rotor, all the blades must be considered simultaneously 
since they will interact, hence a three bladed gimballed rotor would require the nonlinear 
solution of 19,200 to 28,800 DOFs, rendering the analysis prohibitively expensive. 

Additional problems will appear when stability analysis is performed using Floquet’s 
theory, which is standard tool for dealing with the stability of periodic systems. In this 
approach, the stability of the system is assessed from the eigenvalues of the transition 
matrix, which is a fully populated matrix of an order equal to twice the number of spatial 
DOFs. Considering once again the above example, the transition matrix would be of 
order 200 to 300 for a single blade, and 600 to 900 for the gimballed rotor. Furthermore, 
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this matrix is often ill-conditioned because the characteristic frequencies of the system 
vary over an extremely wide range. For instance, if six DOFs are considered at each 
node of the blade, axial frequencies will be included in the model. Such frequencies are 
many orders of magnitude larger than the first flap frequency of the blade, yielding an 
ill-conditioned transition matrix. This limitation is inherent to the approach, and will not 
disappear with increased computational power. 

In view of these numerical difficulties a modal representation appears as a natural 
way of reducing the number of degrees of freedom of the problem. In fact modal 
approaches have been very widely used to analyze rotor blades [3-7,11], and have the 
additional advantage of involving degrees of freedom that have a direct physical meaning. 
However, modal approaches are based on an inherent assumption: the motion of the blade 
is restricted to the superposition of a small number of prescribed modes. Furthermore, 
when applied to nonlinear problems, there is no assurance of convergence or accuracy 
of the procedure. The goal of this research is to develop a finite element based modal 
analysis for rotor blades. The expression finite element based refers to the fact that 
a conventional finite element model of the blade is subjacent to the modal analysis 
which accuracy can be assessed by reference to this complete finite element model. In 
the development of a nonlinear finite element based modal approach, three points are of 
particular importance: the type of nonlinearities, their order, and the choice of the modes. 
The first two point will be addressed in the present chapter and the latter in the chapter 4. 

Consider for instance a nonlinearity of trigonometric type say cosy, appearing in the 
strain energy expression ( 7 is an unknown rotation angle). In the modal approximation, 
this angle is expanded as 7 = Yip', where 7 * are the assumed mode shapes, ip' the 
generalized coordinates, and summation over all assumed modes is implied by the 
repeated indices. To evaluate the strain energy, the expression cos^'ip' must then be 
integrated along the span of the blade; this is of course impossible since ip' are as 
yet unknown, and due to the transcendental nature of the trigonometric functions. To 
avoid this problem, it is customary to expand the cosine function in a truncated series: 
cosYip* « 1 - ^ 7 * 7 iip'xpi . This means of course that the analysis will be limited 
to moderate rotations. Hence, if we wish to develop a modal approximation without 
introducing additional assumptions, the nonlinearities must be of a simple, algebraic type. 

Consider next the order of the nonlinearities, say 7 ", where n is the order of 
nonlinearity. In the modal expansion this becomes YYl k ...ip'ip^ip k ... . It is clear that 
the number of coefficients generated by such expression is proportional to N n , where N 
is the number of assumed modes. Hence, for a 12 mode approximation of an expression 
containing sixth order nonlinearities, 2.9xl0 6 coefficients will be generated, requiring 
22 megabytes of storage on a computer! From this discussion, it is clear that a modal 
approach must be based on expression containing simple, algebraic type of nonlinearities 
only, and the order of the nonlinearities must be as low as possible. 

It is clear that the Reissner’s Principle based formulation described in chapter 2 is 
ideally suited for a modal approximation since it involves nonlinearities only, of quadratic 
order. The details of this modal reduction are given in the following sections. 
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Section 3.2: Modal Approximation to The Total Lagrangian of The Structure. 

When the structure is modeled using the idealizations described in chapter 2 the total 
Lagrangian can be written as: 

C = J U.v: - - (f t c- \f t HF) dV (3.2.1) 

where independent approximations can be made for the displacements, momenta, and 
internal forces. It is convenient to distinguish between rigid and elastic velocities (2.5.21) 
as: 

YX = AV r + K (3-2.2) 

and correspondingly, the following momentum is chosen: 

g = M*AU r +£ (3.2.3) 

Introducing these equations, the total Lagrangian becomes: 

C = C r + C e (3.2.4) 

where the “rigid” Lagrangian is: 

C T = Uj MV_ r ~ \uJ MU r + Uj(C-D) + V t t D (3.2.5) 

with 

M = J A T M* A dV (3.2.6) 

V 

C = j A t M*V^ dV (3.2.7) 

V 

A T jf dV (3.2.8) 

v 

and the “elastic” Lagrangian is: 

c, = J dV (3.2.9) 

The rigid body motions are represented in terms of physical variables 

R t = [Pi P 2 P 3 Qo Qi Q2 Q* J ; R t = [A A A Qo Qi Q2 Qs \ (3.2.10) 

where the P, are the components of the rigid body translation (2.5.7), and Q, the Euler 

Parameters of the rigid body rotation (2.5.1). In contrast, the elastic displacements, 
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momenta, and forces are assumed to be adequately represented by an expression in terms 
of known mode shapes about a given steady reference configuration. 

u, = u? + (3.2.11) 

where u” is the time independent reference configuration of the elastic body, the 
assumed mode shapes, and ifii the generalized coordinates. A similar expansion is 
selected for the elastic forces and momenta, respectively. 

Fi = Ff + Fitf (3.2.12) 

Pi =pf + p\ j 4> i P (3.2.13) 


Section 3.3: Elastic Lagrangian for Beam Elements. 

Introducing the modal expansion described in section 3.2 into the expressions for 
matrices G and H (Appendix A.6) yields: 


G = G l) + G k x!> k u ; H = H° + H k ^ k 


k.k 


where 


and 
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92 
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H° and H k are defined similarly. The elastic rotation matrix becomes 


(3.3.1) 


(3.3.2) 


(3.3.3) 


T = T° + T k j> k u + T kl ^ l u 

where 

T e ° = H^G 07 ; T k = H Q G kT + H k G" 7 ; T kl = H k G ,T 
The strains (2.4.16) become 

§_ = £_ + t£ll>u + e^p k rp[ + e^2jP k xJ} l u ^ 


(3.3.4) 

(3.3.5) 


(3.3.6) 


where: 
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= Tf ai 
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= T kT yH_ 


kirn 


= T klT u' m 


(3.3.7) 
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Similarly the curvature (2.4.16) become: 

K = K° + K k 1p* + 

where 

k 0 = 2G°g'°; K k = 2<?Y* + 2GV 0 ; £ W = 2GV 
the warping strains become: 

e 4 = e 4 + e*t/>* 


Finally the fictitious strains are expanded as: 


eo — £ J £ — 1 — eg + egt4 + CgVuV’ti 


„JW Jfc 


where 


_0 _ OT 0 , t _ « OT t kl _ IT j 

e 9 — 9 • £ - li e 9 — 2? • 9 1 € 9 = 9 ‘9 


The strain expressions (3.3.1) to (3.3.4) can be written in a generic form 


s = £ + M C 


We now turns to the expression of the elastic translation velocity (2.5.10) 

jC = + i:**J + •iViV’l + 


+ + sr* ww 


where 


t£° = T e ° T X° T ^ ; t£* = (7fx* T + T* 7 * 0 ^ ; t>/* = T°V 

= T kT X ,T u£ ; £* w = T ,T u fc ; ; u e * t,m = T ,m ' 


with 


X° = 
and 


X k = 


o - (*0 + U3) (j/0 + U2) 

(20 + U3) 0 -(xo + u^) 

-(yo + U2) (xo + «i) 0 


0 -«3 

0' 


u k 
-«2 


U k 

-u k 


u k 
u 1 


(3.3.8) 

(3.3.9) 

(3.3.10) 

(3.3.11) 

(3.3.12) 

(3.3.13) 

(3.3.14) 

(3.3.15) 

(3.3.16) 

u k 

(3.3.17) 

(3.3.18) 
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The elastic angular velocities are found: 


s£ = tsf + s£ Vi + y'VJ + + si'* Wi 


*k,k 




(3.3.19) 


where 


.0 _ *r(lT. .0 T IT „ 0 . £,•* = 2G » £ t . ^.H = !•« = 2G y 


a;;” = T,”‘ a;” : u> 


The velocity expressions can be combined as: 

yl = vr + j£Vi + + Efy.Vi + 

+ + v^Y'i’ui’iK 

The elastic Lagrangian now becomes: 

£e = £° + £fyj + £}# + £/V>) + £*# + £"^i + ££ #*' 

+ £}U/^ + + Zfd'PW* + c fM^f + C-w'I’Wp 

+ + cOZtfyl*? + £{*>f#tfr 


(3.3.20) 


(3.3.21) 


(3.3.22) 


+ */£:*/* i*:v: + 


klmnikjl j,m:,n 


where 
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(3.3.23) 


Section 3.4: Rigid Lagrangian for Beam Elements 
The matrix M* defined in 2.5.19 can be partitioned as: 

M* = 

The matrix M (3.2.6) can then be written as 

, , i ml 
M = | . _ 

m J 


ml fh* T 

J* 


A r * 

m 


(3.4.1) 


(3.4.2) 
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where 


rh = mX + T e m*T. '} 


J = Xrh T + rhX T - mXX T + T e J*Tj 
The modal expansion of rh gives: 


m = m° + 


where 


^ „ 0 i . ___ k i rpk ^ kl 'T'fc/,’ * 

777 = mu + 1 rn ; m = mu + l e m ; m — i e m 
The modal expansion of J becomes: 

7 = j° + jV* + + J klm ^W^u + J klmn ri<ww 


(3.4.4) 


(3.4.5) 


(3.4.6) 


(3.4.7) 


J° =A' n m or + m°X 0T - mX°X or + T e ° J*T c or 
J* =X°m kT + X*rri 0r + nPX kT + rh k X nT - mX (, X kT - mX k X 0T 
+ J*T^ -|- J* 

J kl =X ( 'm k,T + X k m lT + m k X lT + m kl X (tT - mX k X lT (3.4.8) 

+ T?J*T klT + T kl J*Tg T + T k J*Tl T 
jkim = x k m ,mT + m kl X mT + T k J*T l e mT + T kl J*T™ r 

jklm rpkl j*rpmnT 

The vector Q T — ^ CA T . QE [ T becomes: 

C A = mu + m T wJ + f e m* 

= CA () + CAirpi + + CAl UM + CA$J k u rJ>l 

CB = mu + Ju j? + Xf e m * + T e J*u£ 

= CB ( ' + QBlti + ££}# + C&UWu + QB$JW U 
+ + CBliTu'pWu'Pu'P: + 


(3.4.9) 


CA° = m 0 T y£ ; CA k = m^uf ; CAj = mu k + T*m* = m k 
CA k [ = = (t W + 


if'UY'u 'r t 

(3.4.10) 


(3.4.11) 


CB° = J 0 w? ; ££« = J k ^r ; = mV + X°T k m* + T e °J*^* 

CB kl v = J kl Jl ; CBS = m l u k + X'T*™* + X 1 (r kl + if )m* + T^J*u^ kl + T l e J* u 

CB^™ = = m ,r V + X m (r*' + if )m* + T?J* u* kl + T l e m J*u£ k 

CB kl ™ = J klmn u£ ; CBji = Tl m J*<£ kl 

(3.4.12) 
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Finally the vector D T = | DA?, DB ? J becomes: 

I2A = T e £ 

= DA° + DAM + HAW + DA k JM + DAputfipl + 

(3.4.13) 


and 


where 


DB = X T'£ + 

= DB a + mW, + DBpi’p + + esj i 

+ SbSvMvC + DBXi’ii’ii’Z + 


a4° = T^ 0 ; £ 4 j = T,V" ; = T,V‘ ; £4“ = r t V° i 


and 


DA“ = 

« ; 

n A Mrn 

-■^puu 

+ T?p? ; 

DB k u = 

(x°T k 


r /m *Jfc 
f e E, 


(3.4.14) 


(3.4.15) 




DBt = X 0 T e V * + T e V * ; 


DB 


= ( 


x°r*' + x*rj 


9* 


0 * *0 

+ r e E 


= X l T™£ k 


= ( 


; DBji = (x^ + X l T^£ k + T l e p : 


/„** 

e£i 


X°T e ' m + X 


'37)1 


p‘‘ + X‘"V 
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Section 3.5: Linearization of The Lagrangian 

As a result of the modal expression described in the previous section, the total 
Lagrangian can be written as: 


where 


c = uJmv l - + UJ(C-D) + vJd + c t 
M = M° + M k <f> k + M k, <t> k <t> 1 + M klm 4> k <f> 1 4> m + M klmn <j> k <t> l 4> m <t> n 

c = c Q + c k <f> k + c kl 4 > k <t> 1 + c klm <t> k <f> l <t> m + c* ,m ’vVV n V n 

D = D° + D k 4 > k + D kl <t> k 4> 1 + D klm 4 > k (f> l 4 > m + D klmn <f> k <f, l <f> m <f> n 
C e = C° e + C k <t> k + £ kl <t> k <t> 1 + C\ lm (f> k <t> 1 4> m + C klmn <j> k <f> 1 <j) m 4> n 


(3.5.1) 


(3.5.2) 


The state vector <j> = has been used to simplify the notation. The quasi 

linearization procedure is used to expand the Lagrangian about a known configuration: 

£rr 

£ 6 d 0 

C 




£r 


C = C+ |a£ t , a R r , AU t , A(/> t 

{ 

£r 

£oj 

1 

+ 2 



£4 



+h.o.t 


RR 

£ur 
L C- 4 R 


■'UR 


Cuu 


£<t>R ^ 4 >U 



AR 


AR 


AU 


A <f> 


} 


(3.5.3) 



where the first derivatives of the Lagrangian are: 


£r = Z(MU r + D) ; C k = Z{MTLr + D ) ; 

Cu = (MV r + C)~ ( MU t + D) ( 3 . 5 . 4 ) 

£4, = UjM t V r - | UjM t U T + +(]£?- 

The second derivatives of the Lagrangian are: 

£rr = W ; £ RR = W ; £ur = MZ t ; £ UR = MZ T ; £jjjj = -M 
£<t>R = (vjMi+Df)z T ; £ 4>k = (ujM,+D[)z T ] 

C<t>u = (MZr + £,) - (MiU r + Bi ) ; (3 - 5 - 5) 

tn = - l -UjM tJ U r + + (y-r - Uj)D tJ + £ e ij 

Increments in the rigid velocities can be related to increments in the rigid body parameters 
as: 


AV r = Z l AR + Z T AR 
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The following quantities are now defined: 

C k = C k + Uj M k V_ r - i UjM k U r + Ujc k + (vj - Uj)D k 
C kl = C kl + Uj M kl V_ r - i UjM kl U r + UjC kI + (vj - Uj'j D kl 
£klm = £ Wm + uJ M klmy r _ hjJ M klm U_ r + Uj C klm + (vj - Uj^jL 

(vj-l 


c klmn = c klmn + uJm 


klmn V T - - l -UjM klmn U T + lljc 
& 


T / ~<klmn 


+ 


and 


- o 


^klm 


c =M"v r + Q.; 

& = M k V r + C k ] 
C kl = M kl V r + C kl -, 


if = M n u r + D° 
if = M k U r + D k 
D U = M kl U r + D kl 


b ktm = M klm U T + D klm 


b kimn = M klmn 


Ur + D 


i klmn 


cr m = M kim v r + c um ] 

jjkhnn ^ ^jklmny _j_ gklmn . 

The first derivative becomes: 

Cr = ZY; C k = ZY ]Cu = 2 L-Z; C 4>i = C et 
The second derivative becomes: 

Crr = W ; Crr — W ; Cur — MZ ^ = MZ ^ ; Cuu ~ —M 

C<t>R — Y.{Z t ; C^r = ¥. t Z T ; C<t>u = - Z, 

C<j><t> = C-eij 

where 

x = a ° + c V + £‘W + 

r = i>° + + £‘Vy + £“”VVV m + fi"" v*Vv 

The following arrays were defined 
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Chapter 4 

Comparisons of Modal and Finite Element Methods 
Section 4.1 Introduction 

The appropriate choice of modes is crucial to achieve accuracy in the modal analysis. 
In general, natural vibration modes have been selected in modal analysis of rotor blades. 
The relative merits of various sets of modes have been investigated, for instance, coupled 
or uncoupled free vibration modes of a rotating or non-rotating blade [4-1,3]. It 
is important to note that natural vibration modes characterize the linearized dynamic 
behavior of the blade, i.e. the dynamic behavior of small, time dependent perturbations 
about a given, steady equilibrium position of the blade. Even though it is natural to use 
such modes in the analysis of nonlinear problem, it is well known that the accuracy and 
efficiency of a modal method depends on the “quality” of the assumed modes, i.e. the 
ability of the assumed modes to represent the actual response of the blade. 

When natural vibration modes are used in conjunction with a displacement based 
energy formulation that includes axial displacement as an independent variable, the 
performance of the modal approximation is extremely poor. Consider the lateral deflection 
of a blade in the nonlinear range, under a simple tip oscillatory load. If flapping modes 
are used in the modal approximation, the lateral deflection is found to be much smaller 
than that predicted by the full finite element model. This can be explained by the fact 
that flapping modes contain no axial component (since they are linearized modes), hence 
foreshortening of the blade is not allowed in the modal approximation and this results in 
large axial loads which in turn, stiffen the blade considerably. The situation is somewhat 
improved by adding axial vibration modes, but a large number of these modes is required 
to obtained a good solution. 

The reason for this behavior is twofold. First, in a displacement based formulation, 
the stress-strain relationships are strongly enforced (i.e. they are satisfied on a point by 
point basis), therefore, a very small error in the estimation of the axial strain (as should 
be expected from a modal approximation) will result in very large axial forces, because 
of the very large axial stiffness of the blade. In fact, the inextensibility assumption is 
often made to avoid this problem, however, the formulation is then restricted to single 
load path blades. This problem can be overcome when using the mixed formulation 
described in Chapter 2. Indeed, in a mixed formulation, the stress-strain relationships are 
only enforced in a weak sense (i.e. in an integral sense); hence, small errors in strain do 
not necessarily result in large errors in the forces. 

Second, the actual axial displacement of the blade is due to foreshortening (a 
nonlinear kinematic phenomena), whereas axial vibration modes characterize true axial 
vibrations (a purely linear vibratory phenomenon). In other words, we are trying to 
"synthesize" a nonlinear kinematic mode shape, with linear vibratory mode. Since these 
two phenomena are not physically related, we should hardly expect to obtain good results 
in predicting axial displacements. This discussion has focussed on axial displacements 
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due to foreshortening however, the above arguments equally apply to any nonlinear 
behavior of the blade. For instance, transverse loads applied to the blade create a torque 
due to the blade’s transverse deflections. This nonlinear kinematic coupling is very 
important for helicopter blade response, as it can change its angle of attack. 

This clearly indicates the need for selecting alternate mode shapes that contain 
information about the nonlinear behavior of the structure. Several concepts have been 
proposed to improve the quality of the modal bases when dealing with nonlinear problems. 
The conceptually simplest method it to recalculate a new set of natural vibration modes 
every once in a while as the deformations of the blade become significant [4-7], In fact, 
the natural vibration frequencies and associated mode shapes of a helicopter blade are 
known to vary significantly around the azimuth [4-3]. Even though this approach might 
give good results, it does so at a tremendous computational cost, since the modal basis 
must be updated during the response calculation, and the modal reduction scheme must 
be repeated at each update. Another approach is to include in the modal basis natural 
vibration modes about various different equilibrium configurations of the structure. This 
method is attractive since only a modest increase in computation cost is required to 
evaluate the various equilibrium configurations. Furthermore, this method appears to 
give accurate results, see for instance [4-7], 

Finally, the concept of perturbation modes seems to hold promise for improving the 
accuracy of modal methods. It was introduced by Thompson and Walker [4-4] as an 
analytical tool for the study of the nonlinear behavior of beam structures, and later the 
same concept was used by Noor et al [4—5,6] in the nonlinear static analysis of beam and 
shell structures in conjunction with the finite element method. The very same concept 
will be used here to study nonlinear dynamic problems. 

Section 4.2 Perturbation Modes. 

Static perturbation modes can be evaluated from a finite element model according to 
the following procedure. The incremental form of the finite element equation is: 

K(u)Au = Q - F(u), (4.2.1) 

where u is the vector of nodal unknown that includes both nodal displacements and 
forces, K the stiffness matrix linearized about a reference configuration u, Q the vector 
of externally applied loads, and F the vector of equivalent nodal forces. Equilibrium is 
achieved when A u = 0 , or 


F(u*)-Q = 0, (4.2.2) 

which simply states that at equilibrium, the equivalent nodal forces are equal to the 
applied loads. The equilibrium configurations* is of course a function of the applied 
load. Consider now an applied load of the form A Q (A is a scalar), the equilibrium 
condition is: 


F t (u(\))-XQ, = 0. 


(4.2.3) 
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Since F, is a nonlinear function of the displacements, equation (4.2.3) can not be solved 
easily, however, taking the first derivative of this relation with respect to A yields: 


dFj 

duj 



= Qu 


(4.2.4) 


where = K tJ is the linearized stiffness matrix, andu^ is the first perturbation mode 
which is clearly nothing else but the solution of the linearized problem. The second and 
higher order derivatives become: 


i j U 



d 2 F, (i) (i) 
dujdu k 3 ’ 


(4.2.5) 
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du,du t ( 4u > 3> “‘ 1 + 3 “f )u i 2> ) 
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d*F, 


dujdu k dm 


U - U k Uf , 
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.^L( 5 ^'u»> + 10 A< 21 ) 

dujduk \ 3 k 3 k ) 

^ ri0u (3) u (1) u (1) + 15 u (2) u (2) u (1)> 1 

dujdu k du t \ w 3 k 1 + 5 3 k u ' r 


d*F t 


(4.2.6) 


(4.2.7) 


(4.2.8) 


These relationships are recursive, and involve a single inversion of the linearized stiffness 
matrix. They also involve higher order derivatives of the equivalent load vector F t . This 
task is relatively simple when dealing with the finite element formulation described in 
Section 2, since the energy expressions are purely algebraic, quartic expressions. This 
also explains why fourth and higher order derivatives vanish and arc thus absent in (4.2.7) 
and (4.2.8). In a perturbation theory approach, the solution would be written as: 


«, = Hi + A*! 1 ’ + ^A‘ u ;" + ^A'X"' + 


»..(»> 


k 3 < 3 ) 


2! 


3! 


(4.2.9) 


however, the convergence characteristics of this expansion arc extremely poor. A much 
better approach is that proposed by Noor et al. [4-5,6] where the perturbation modes are 
simply added to the modal basis of a standard modal analysis as described in Chapter 3. 

The above formulation is limited to static problems; however, it can be readily 
extended to accommodate dynamic situations. Let Q be the inertia forces associated with 
a natural vibration mode shape u,, i.e. Q = where M is the mass matrix, and u 

the associated natural frequency. The recursive relations (4.2.4) to (4.2.8) can be used to 
obtain perturbations of these natural vibration mode shapes. Such modes will be termed 
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here “dynamic perturbation modes.” In this case the stiffness matrix is the sum of the 
structural stiffness matrix, and the centrifugal stiffness matrix. 

Section 43 Numerical results and discussion. 

The purpose of this chapter is to assess the accuracy of modal analysis methodologies. 
This will be done by computing the dynamic response of structures obtained on one hand 
from the full finite element procedure, which will be taken as a “reference” solution, and 
on the other hand, from the modal analysis procedure , with various modal bases. The 
computation proceeds with the following steps. First, the physical structure is discretized 
into a number of beam elements and the corresponding finite element equations are 
integrated in time using the finite element in time procedure to obtain the reference 
solution. The second step is the selection of a modal basis consisting of a mixture 
of the following types of modes: natural vibration mode shapes about the reference 
configuration, natural vibration mode shapes about any other configuration, and static or 
dynamic perturbation mode shapes. The third step consists of the modal reduction. It is 
important to note that the full finite element model, the modal basis, and the modal 
reduction are all based on the identical finite element discretization of the physical 
problem. In the last step, the modal equations are integrated in time to obtain the modal 
response. In all cases, both full finite element and modal equations are integrated using 
two noded elements in time (i.e. a linear approximation for the displacements within 
each time step), and identical time step size are selected. 

It is important to note that all the models discussed here are based on the exact same 
equations, namely the Euler equations resulting from the minimization of the Lagrangian 
expression. The only difference among the various solutions is the description of the 
solution fields: in the full finite element model, the solution is represented by polynomial 
expressions defined within each finite element, whereas in the modal analysis, the solution 
is represented by the modal superposition. Hence, all the responses presented in this work 
are based on the exact same equations, with different spatial discretization of the solution. 

The first test case consists of a straight, cantilevered blade, with a thin-walled, 
rectangular cross-section. The blade has a length of 3 m, a width of 0.15 m, and a 
height of 0.02 m. The wall thickness is 1mm, and the material is aluminum (Young’s 
Modulus 73 GPa, density 2700 kg/m % ). The overall geometry of the blade is depicted in 
Figure 4.2. The blade does not rotate, and is subjected to a tip load of 250 N oscillating 
with a period of 1 second. The blade is modelled with four cubic beam elements, for 
a total of 96 displacement degrees of freedom. Forty time steps are used to model the 
1 second period. Table 4.1 details the three different modal bases used for correlation. 
Figure 4.3 compares the flapping deflections of the tip of the beam for the various modal 
bases and the full finite element model. A good correlation is obtained for all modal 
bases, even though very large transverse deflections occur (1.3 m compared to the 3m 
length of the beam). Figure 4.4 shows the correlation for the axial displacement, i.e. 
the foreshortening of the tip of the beam. Note that a single perturbation about the first 
flap mode yields an excellent correlation, whereas adding the first 3 or 5 axial vibration 
modes do not achieve this level of accuracy (bases 2 and 3, respectively). 
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Table 4.1 Description of the Modal Bases for the 0° Nonrotating Blade 



Basis 1. 

Basis 2. 

Basis 3. 

Flap Modes 

1 

1 

1 

Flap Perturbations 

1 

0 

0 

Lead-Lag Modes 

1 

1 

1 

Axial Modes 

0 

3 

5 


In the second test case, the blade’s cross section is tilted at a 45° angle with respect 
to the loading axis. Table 4. 2 summarizes the various modal bases used here for the 
modal analysis. Note that basis 4 involves modes which were taken about the deformed 
configuration of the blade under the static load of 250 N at the tip. Figures 4.5 and 4.6 
show the in-plane and out-of-plane deflections of the blade, which are all in reasonable 
agreement with the reference solution. Figure 4.7 shows the tip twist of the blade, a 
nonlinear kinematic phenomenon due to the offset of the tip load creating a torsion 
moment arm. Basis 4, which contains the natural vibration modes about a predeformed 
configuration of the blade performs well when the dynamic response of the blade is in 
the same direction as that of the predeformation (the first half of the period), however 
it performs very poorly when the dynamic response is in the opposite direction of the 
predeformation (the second half of the period). This clearly shows that modes about 
a predeformed configuration should be avoided when the dynamic response involves 
complete reversals, as is the case for a helicopter blade. Basis 3 contains natural vibration 
mode shapes only, and performs very poorly, missing the tip elastic twist by over a factor 
of two, even though 5 torsion modes were used in an attempt to capture this kinematic 
phenomenon. The reason for this poor correlation is that the observed twisting of the 
blade is due to a nonlinear coupling effect, whereas the natural vibration mode shapes 
characterize true torsion vibrations. These two phenomena are not physically related, 
and this explains the poor correlation. Finally, Figure 4.8 depicts the amplitude of the 
torsional warping deformation. Only basis 1 provides a good correlation for this quantity 
that is directly related to the torsional loading in the blade. Clearly the perturbation 
modes of basis 1 outperform all other bases, even though it only includes a total of 5 
modes. Note that a static perturbation mode was included in this basis to provide the 
proper nonlinear coupling between transverse loading and twisting. 
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Table 4.2 Description of the Modal Bases for the 45° Nonrotating Blade 



Basis 1. 

Basis 2. 

Basis 3. 

Basis 4. 

Flap Modes 

1 

1 

1 

1 

Flap Perturbations 

1 

1 

0 

1 

Lead-Lag Modes 

1 

1 

1 

1 

Axial Modes 

0 

0 

5 

0 

Torsion Modes 

0 

5 

5 

0 

Static Modes 

1 

0 

0 

0 

Static Perturbations 

1 

0 

0 

0 


In the third test case, the 45° blade is now spinning at an angular velocity of 6.28 
rad/sec, and is subjected to a 350 N oscillating tip load. Table 4.3 summarizes the various 
modal bases used in this case. Figures 4.9 and 4.10 show the in-plane and out-of-plane 
deflections, which are all in good agreement with the reference solution. Figures 4.11 
and 4.12 show the tip twist and torsional warping amplitudes. Once again, basis 1, which 
involves perturbation modes, clearly outperforms the other bases, even though it includes 
5 modes only. 

Table 4.3 Description of the Modal bases for the 45° Rotating Blade 



Basis 1. 

Basis 2. 

Basis 3. 

Flap Modes 

1 

1 

1 

Flap Perturbations 

1 

1 

0 

Lead-Lag Modes 

1 

1 

1 

Axial Modes 

0 

0 

5 

Torsion Modes 

0 

5 

5 

Static Modes 

1 

0 

0 

Static Perturbations 

1 

0 

0 


The last test case involves an actual helicopter blade: Sikorsky Aircraft’s Blackhawk 
blade. This 27 ft long blade is modelled with 16 cubic elements, for a total of 336 
degrees of freedom. Forty time steps are used to model a single period of 0.23 seconds. 
The aerodynamic loading is approximated by a concentrated lift (1000 lb) and drag 
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(250 lb) forces applied at 92% span. Figures 4.13 and 4.14 show the in-plane and out- 
of-plane deflections of the blade, which are in excellent agreement with the reference 
solution. Figure 4.15 shows the tip twist (i.e., the tip angle of attack) of the blade, and 
large discrepancies are observed between the reference solution and the various modal 
responses. Several bases are examined, but produce only marginal improvement. The 
reason for this discrepancy is probably the presence of nonlinear coupling between the 
twisting of the blade and rotational dynamic effects. Such nonlinear couplings are non 
properly represented by natural vibration modes, nor by perturbation modes. Indeed, 
when calculating the perturbation modes, gyroscopic terms are ignored. 

It is important to note that all the test cases examined in this effort involve a prescribed 
loading. In actual problems, the loading is of an aerodynamic origin, and hence dependant 
on the response of the blade, most noticeably on the angle of attack. Were the above 
modal analysis used in an actual coupled problem (aerodynamics coupled with structural 
dynamics), the discrepancy observed in the angle of attack (Figure 4.15) would generate 
different loading conditions, which in turn would further change the blade’s response. 
This would generate different responses for flapping, lead-lag, and twisting. 

Table 4.4 Description of the Modal bases for the Blackhawk Blade 



Basis 1. 

Basis 2. 

Basis 3. 

Flap Modes 

3 

3 

3 

Flap Perturbations 

0 

0 

1 

Lead-Lag Modes 

2 

2 

2 

Axial Modes 

0 

1 

1 

Torsion Modes 

5 

5 

5 

Static Modes 

1 

0 

0 

Static Perturbations 

1 

0 

0 


Finally, it is interesting to compare the computational times for the various ap- 
proaches. Table 4.5 summarizes the CPU times for the full finite element analysis and 
the various modal approaches, normalized by the CPU time for the full finite element 
analysis. It is interesting to note that even though the 5 mode modal analysis only re- 
quires a small fraction of the full finite element CPU time, the cost of the modal analysis 
drastically increases with the number of modes. In fact, the 12 mode analysis is more 
expensive than the full finite element model. As the complexity of the full finite element 
model increases, its cost will increase as well, however. Table 4.5 clearly indicates that 
the costs of dealing with modal or full finite element models become comparable as the 
number of modes increases. Since the accuracy of the modal analysis is questionable even 
when using an increasing number of modes, the full finite element model is preferable. 
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Table 4,5 Normalized CPU times for the 45° rotating blade 


Analysis Method 

Normalized CPU Time 

Full FEM (96 DOFs) 

1.00 

Basis 1. (5 Modes) 

0.090 

Basis 2. (8 Modes) 

0.261 

Basis 3. (12 Modes) 

1.22 
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Fig 4.3. Axial Tip Displacement vs. Time for the Nonrotating 0° Blade. 
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Fig 4.4. In-Plane Tip Displacement vs. Time for the Nonrotating 45° Blade. 
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Fig 4.6. Tip Twist vs. Time for the Nonrotating 45° Blade 
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7. Tip Torsional Warping vs. Time for the Nonrotating 45° Blade. 
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Fig 4.8. In-Plane Tip Displacement vs. Time for the Rotating 45° Blade 
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Fig 4.9. Out-of-Plane Tip Displacement vs. Tune for the Rotating 45° 
Blade. 
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Fig 4.10. Tip Twist vs. Time for the Rotating 45° Blade. 










Fig 4.11. Tip Torsional Warping vs. Time for the Rotating 45° Blade. 
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Fig 4.12. In-Plane Tip Displacement vs. Time for the Blackhawk Blade. 
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Fig 4.13. Out-of-Plane Tip Displacement vs. Time for the Blackhawk Blade. 


O 



Fig 4.14. Tip Twist vs. Time for the Blackhawk Blade. 
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Chapter 5 


Kinematic Constraints between Elastic Bodies. 

Section 5.1: Introduction 

This chapter deals with the formulation of the kinematic constraints between the 
elastic bodies of a multibody system. 



To define the multibody system, the initial position and orientation of the various 
elastic bodies must be given, i.e. p('\t = 0) the inertial position of the material point 
O (2.5.7), and Tr'\t = 0) the inertial orientation of the body (2.5.1) are specified. It is 


5 - l 



assumed that T r °^(i = 0) is an identity matrix. The notation (.)^ is used to identify the 
body i whenever the distinction is necessary. 


Section 5.2: Displacement constraint hinge 


A hinge corresponds to a set of kinematic constraints on the relative displacements 
and rotations of two distinct bodies, and involves the degrees of freedom of two nodes, 
one in each body. A set of orientation vectors d, is also associated with each hinge to 
allow for the definition of the relative rotation constraints. 

Consider first the kinematic constraint corresponding to the continuity of displace- 
ments across the “hinge” which can be written as: 


r[ 1} = R? } 


(5.2.1) 


where R' h is the position of the hinge calculated in body i. In this section we examine 
the formulation on the problem when constraint 5.2.1 is enforced. However, it is also 
possible to enforce the time derivative of this constraint, this approach will be discussed 
in section 5.4. When the constraint conditions are adequately formulated, the Lagrange 
multipliers become the unknown forces transmitted at the “hinge” Equation (5.2.1) can 
be expanded as: 


(P {1> + 4h + 4 1} ) - (^ (2> + + 4 2> ) = 0 (5.2.2) 

where P is the position of the reference frame (see (2.5.3), r 0 j, the position of the 
hinge in the undeformed body, and the elastic displacement of the hinge. The vector 
relationship (5.2.2) can be expanded in component form as: 


Ui h h J 
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(5.2.3) 


The reference triad of each elastic body is related to the inertial triad through (2.5.2), 
hence the constraint condition becomes: 


h,h,h\{ 
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or, in component terms: 

p(l) + T T {1} Tr {1) u[ l) - pW _ y r <2> Jjp(2) (2) _ 0 


}=0 

(5.2.4) 

(5.2.5) 
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Consider now the kinematic constraints on the relative rotations of the two bodies. 
At time t = 0, the orientation of the hinge is given as: 


50) 


d 


h 

% 


d 



(5.2.6) 


This orientation triad can be related to the base vectors in each body using (2.5.2) and 
(2.2.7): 
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At time t, the two bodies have now rotated with respect to each other: d^ becomes 
d*^, however, d*^ is attached to body 1, and d*^ is attached to body 2, hence: 
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With the help of (2.2.14) and (2.5.2), we find: 


d?"> 

= t (a)T T (a)T T 0{ a )T T (a)T 
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= r^ t 
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(5.2.8) 


(5.2.9) 


Let d\ a \ gj a \ Q^ a \ q{°\ and be the Euler Parameters associated with the 
rotation matrices T$ a \ Tr^ a \ T r , and R respectively. The relative rotation 
rule (A36-37) yields: 

r<"> = (5.2.10) 


The relative rotation at the hinge is now: 


dT 

4 (1> 

= rWrW 

d'< 2 > 

dl< 2 > 

= 5 

£< 2 > 

£< 2 > 
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df> 


4< 2 > 


and the relative rotation rule (A38) now yields the Euler parameters associated with S 
as: 

Si = r< 1)T SV< 2 > (5.2.12) 
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where 


S 1 


s 2 


s 


3 


0 10 0 

-10 0 0 
0 0 0 -1 
0 0 10 

0 0 10 

0 0 0 1 

-10 0 0 
0-100 

0 0 0 1 
0 0-10 
0 10 0 

-10 0 0 


(5.2.13) 


(5.2.14) 


(5.2.15) 


The appropriate relative rotation can be constrained by enforcing the corresponding Euler 
Parameter to vanish, i.e. 


Si = 0 


(5.2.16) 


The kinematic constraints (5.2.5) or (5.2.16) can be enforced by means of the 
Lagrange multiplier technique. Since a Hellinger-Reissner formulation is used for the 
other components of the Lagrangian, it is convenient to use a similar formulation for 
the constraints, namely: 


h = - p< 2 > - r r < 2 >r«uf>) 

+ rJ 1 > T Sr_m-l.(f. L + gT s ) 
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and gi are the Lagrange multipliers, and a is a large number. 


(5.2.17) 


(5.2,18) 


Section 5.3: Modal approximation of the kinematic constraint 
at a displacement hinge. 


In this section the kinematic constraint at a hinge (5.2.17) is expanded using the 
modal approximation described in section 3.2. The elastic displacement of the hinge 
(5.2.5) is expanded as: 


and 
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where 


and finally: 


where 
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Similarly, the array r( a ) (5.2.10) is expanded as: 

r<*> = r°<“> + 

where * 

r°< a > = 

L k (o) _ ^4 A B (d^q k < a \ 

With the help of (5.2.1) to (5.2.7) the constraint condition (5.2.17) becomes: 

H = f (£<>> - P< 2 > +£*(»- E *») + r - -L (/./ + /.p) (5.3.8) 

The constraint expression is a nonlinear funtions of the rigid body parameters and 
generalized elastic coordinates of the two bodies. The quasilinearization procedure can 
be used to expand this expression about a known configuration to find: 
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The derivative with respect to is: 
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where the following notation is used: 

sign(l) = +1; sign{2) = —1 
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(5.3.12) 

(5.3.13) 
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The derivative with respect to Vv is: 
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The derivative with respect to BP\R^ is: 
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The derivative with respect to is: 
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Finally, the derivative with respect to /, / is: 
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Section 5.4: Velocity constraint hinge 

The second constraint is the continuity of velocities across the “hinge” which can 
be written as: 

U<1> U<2> 

Rh = Rh (5.4.1) 

where is the velocity of the hinge calculated in body i. When the constraint 
conditions are adequately formulated, the Lagrange multipliers become the unknown 
impulses transmitted at the “hinges”. Equation (5.4.1) can be expanded as: 

'pm + 4> + 4 l >) - (?\» + 4> + af >) = 0 

The vector relationship (5.4.2) can be expanded in component form as: 
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The reference triad of each elastic body is related to the inertial triad through (2.5.2), 
hence the constraint condition becomes: 
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The kinematic constraints (5.4.4) or (5.4.5) can be enforced by means of the Lagrange 
multiplier technique. Since a Hellinger-Reissner formulation is used for the other 
components of the Lagrangian, it is convenient to use a similar formulation for the 
constraints, namely: 
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Section 5.5: Modal Approximation of the kinematic constraint at a velocity hinge. 

In this section the kinematic constraint at a hinge (5.4.6) is expanded using the modal 
approximation described in section 3.2. The elastic displacement of the hinge (5.4.5) is 
expanded as: 
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and 
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With the help of (5.5.1) to (5.5.5) the constraint condition (5.4.6) becomes: 

h = f (£<»> - pm + „;<*> _ „;< 2 > + „•<» _ „f > + „;<■> _ „;<’>) _ _L If A 

2 “ V (5.5.6) 

The constraint expression is a nonlinear funtions of the rigid body parameters and 
generalized elastic coordinates of the two bodies. The quasilinearization procedure can 
be used to expand this expression about a known configuration to find: 
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The derivative with respect to R ^ is: 
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The following notation is used: 


sign(l) = +1; sign( 2) = -1 
(0) = 2 when (a) = 1; (0) — 1 when (a) = 2 


( 5 . 5 . 15 ) 


( 5 . 5 . 16 ) 


The derivative with respect to is: 

Ku'M = sign{a)fj + v*^') 


( 5 . 5 . 17 ) 


The derivative with respect to is: 

K u ,\ a) = sign{a) (/ r u" <o) ) 


( 5 . 5 . 18 ) 


The derivative with respect to / is: 


Hf = 


P <1> -P <2> +v0 


+<1> - - v$ <i> ^ ^ <2> - -*<'> 


+ ”3 -113 -//a 

( 5 . 5 . 19 ) 


The derivative with respect to R^ Q \R^ is: 


H R ( a ) R { a ) = sign(a) 
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The derivative with respect to R^°\R^ is: 


Hrmjkc) = sign(a) 
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The derivative with respect to ipP, R^ is: 
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The derivative with respect to rpu , is: 


rp 

H u i {a)R{a) 


sign{a) 


0 

0 
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hc\ M + /2C'<“> + / 3 cf> 
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The derivative with respect to f,R^ is: 
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The derivative with respect to f,R ^ is: 

Hf R ( a ) = sign(a) 
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The derivative with respect to /, ipu*'* is: 


H fui{a) = sign(a) 

The derivative with respect to is: 

H fu , ia) = sign (a) 
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( 5 . 5 . 30 ) 
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( 5 . 5 . 32 ) 
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Finally, the derivative with respect to f,f is: 


H ff = - 


1_ 

a 


I 0 

I 


(5.5.33) 


Section 5.6: Finite Element in Time Discretization of a Single Elastic Body. 


The finite element in time procedure will be used to derive the governing equations 
of the problem. In the modal formulation, the time varying unknowns of the problem 
are the dispacement variables U_, the force variables £, the momenta variables V, and A. 
The displacement variables are interpolated as follows: 

N u 

Jt=l 


where g * are the interpolation functions, and U^_ the nodal displacements. If t x and if 
are the initial and final times of the period under investigation, the nondimensional time 
t is given as: 



where A t — tj — t x . The shape function will be chosen as follows: 

si = \(Po - Pi) 

si = j(Po + Pi) (5.6.3) 

si = Pt ^ k ~_ P ^ = Si- 2 i = 3,4,...W, 


where the P, are the Legendre Polynomials, and the S t their integrals. In view of the 
properties of Legendre polynomials, these interpolation functions are such that: 


*1(-1) = 1.0; gl(- 1) = 0.0; g k u (- 1) = 0.0, k = 3, ...N u 
*1(+1) = 0-0; gl(+l) = 1.0; ff J(+l) = 0.0, k = 3 ,...N U 


which means that: 


UiU) = ^i,andW(f/) =U 2 


(5.6.5) 


The time derivative of the displacement generalized coordinates are now: 


Nu 


jfc=i 


(5.6.6) 
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where 



•2 1 
<7u = 2 

9u = Pk-2 k = 3,4,...N u 


(5.6.7) 


The force and momentum variables arc interpolated with the help of the Legendre 
polynomials as well: 

N, 

£(*) = (5.6.8) 

Jfc=i 

and 



a<) = E^( r )& 
*= 1 

(5.6.9) 

where 


g h f = P k -i A: = 1,2, ...N f 

(5.6.10) 

and 


ffp = Pk-i k = 1,2, ...N p 

(5.6.11) 


These various interpolation functions are conveniently written in matrix form, for 
instance (5.6.1) can be written as: 


U{t) = B u ti- U{t) = B V U;, 
£(t) = B f T ; V{t) = B P V 


where 


ul = [af,& T al\ 


(5.6.12) 


(5.6.13) 


The complete Lagrangian for a single elastic body was derived in chapter 3, equation 
(3.5.3), introducing the interpolation functions we find: 


£ = £+[ AU T A j? AV T J { 


1 

+ 2 


£uu Out 

Cut 

AU 

£tt 

0 

At 


C-pv 

AT 


4* 

Lt 

Lv 


(5.6.14) 


} + h.o.t 
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16 


where 


and 


</ 

4 = j (bZCu + BSCv)*; 

t, 

if 

jrdt- (5.6.15) 

t, 

if 

■pdt 




ij 

Buu = j (b^CuuBu + B T U C UV B V + B^Bu V B u ^jdt\ 

t. 

l s 

But = J B^CuTBfdt ; 

t, 

if 

But = J (b^ But B p + B„ BvvBpjdt-, 

t. 


if 

Btt — J BjCjrjrBfdt ; 

u 


Btt = 


if 

J BpCTTBpdt 

t. 


(5.6.16) 


In the case of mutually interacting bodies, kinematic constraints such as (5.2.2) 
will couple the behavior of the various bodies. However, these constraints are purely 
kinematic, i.e. they do not involve forces and momenta, hence, forces and momenta are 
independent variables within each body. The stationarity condition with respect to these 
variables yields: 


A £ = -CZ^jrBjr - CttBujtAU 


'TT±T 


-1 rT 


tt^ut’-^- 


A — — BjJpB'p — C^Cff-pAU 

Introducing these results into (5.6.14) yields a reduced Lagrangian expression: 

B = B + AtfLu + \a^L U uAU + h.o.t 

£ 


(5.6.17) 


where 

and 


\=u — By — ButB^Bj: — ButB"p\>B-p 
L uu = Buu - Bu?Bf l jrBuj: - BuTBp\ > Bj iV 


(5.6.18) 

(5.6.19) 

(5.6.20) 
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Section 5.7: Finite Element in Time Discretization of a Hinge Constraint. 

Consider a single hinge joining two elastic bodies. The corresponding kinematic 
constraint is in the form of (5.2.1). The Finite Element in Time discretization of 
the displacement variables was given in section 5.6, and the Lagrange multipliers are 
interpolated as: 


N u 


/£(<) = 

Jt=i 


(5.7.1) 


where 

9i(r) = Pk-l A: = 1 , ...Np (5.7.2) 

and the corresponding interpolation matrix is B^. 

With the help of these interpolation functions, the constraint condition (5.2.9) be- 
comes: 

H = H + [ A U^>' r A U^l_ A/) 7 ’] { 'Hym 


n u 


i 

+ 2 


^uo)u< o Wi/wum n U ( i)^ 

Am 

'Hitwuw "Hymn 

A 7/< 2 > } 

n^i 

A£ 


(5.7.3) 


} + h.o.t. 


where 


l f l / 

n m = J 7 % = J Bln t 


dt 


(5.7.4) 


and 


t, 

tj 


'tiuwuoy = j Bl B u dt\ n U (,) M = J Bln^^Bfidt; 

U 

if 

n w = J Bln^B^dt 


(5.7.5) 


The Lagrange multipliers are independent variables for each hinge. The stationarity 
condition with respect to these variables yields: 

= -Kin, - - n^n^Au^ (5.7.6) 

Introducing these results into (5.7.3) yields a reduced constraint expression: 


H = H + [ A AU^E\{ 




H W { 2 > 


i 

+ 2 




a 

A UW 


(5.7.7) 


} + h.o.t 
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where 


(5.7.8) 


and 




(5.7.9) 


Section 5.8: Governing Equations of Multibody Systems. 


Consider now a problem involving jVj elastic bodies interacting through Nh hinges. 
Hamilton’s Principle writes: 

6L + 6W + 6{^[U(ti)-Vi\ -P£[u(*y)-Uy]} =0 (5.8.1) 

where L is the total Lagrangian of the system given by: 

N b Nh 

l = ^2c {t) + j2 H{a) ( 5 - 8 - 2 ) 

1=1 0=1 

8\N is the virtual work done by the applied loads; F\ and Py are the initial and final 
momenta vectors, respectively; Uy and Uy the initial and finaTvalues of the displacement 
vector U(<); and finally t t and fy the initial and final times of the period under 
consideration. The unknowns of the problem are the displacements in each of the jVj 
elastic bodies, i.e. 





(5.8.3) 


With the help of the finite element in time discretization described in the previous sections, 
the unknown nodal values are: 


fiT _ 1 T T f i T 1 

[ U 1 i u 2^— 

(5.8.4) 

0?-= i = l,2 

(5.8.5) 


where 


The variations of (5.8.1) with respect to the unknown initial and final momenta yields: 


Oi=Oi ;U2 = 0y (5.8.6) 

and the variation with respect to 0 yields: 

U(t t ) = U i andU(< / ) =Uy (5.8.7) 


which, in view of (5.6.5) is simply: 


The variation of (5.8.1) with respect to the unknown displacements yields: 


F + KAU+P +Q= 0 


(5.8.9) 


where F is the equivalent load vector and K the stiffness matrix, both obtained from the 
appropriate assembly of the Lagrangian (5.6.20) of each elastic body and the constraint 
(5.7.7) of each hinge. The momentum vector is: 



pr,-pJ,o,...o| 


(5.8.10) 


and Q is the vector of applied loads. 

The intermediate nodal variables can be evaluated and eliminated from (5.8.9) to 
give: 


K„ K if 

AUy 


-Pi 


F. + Q, 

K // 

AU/ 


+P / 


F/ + Q/ 


In a step by step integration scheme the initial displacement and momentum are known 
quantities, hence AU^ = 0 and the corresponding quantities at the end of the time step 
are: 


ALV = -K- / 1 [P 1 + F i +Qj 
P/ = Ky/AU/ + F / + Q / 


Section 5.9: Numerical results and discussions 
Section 5.9.1 Spinning Top 

The spinning top will be investigated to test the modeling of rigid body motions 
and constraint equations. The base point of the spinning top is not free to translate 
but rotations are allowed. As discussed in the previous sections, this constraint can be 
applied in two ways: the displacement of the base point can be constrained ( we will 
refer to this case as the “D-hinge”.), or the velocity of the base point can be constrained 
( the “V-hinge”). The numerical accuracy and stability of these two models will be 
addressed here. 

The analytical solution of the spinning top case is obtained by using Euler angles 
6, <f>, and ip. The Lagrangian of the system writes: 

L = T - U 

= ~ PQdcosd (5.9.1) 

where /„ and u t are the moments of inertia and angular velocities in 1,2 and 3 
coordinates respectively. ji,g, and d represents mass density, gravitational acceleration 
and the distance from the origin to the center of gravity of the body respectively. By 
using the Euler angles 8, <p, and ip Eq. (5.9.1) can be rewritten as: 

L= -Ii<^ 2 sin 2 0 + -I 2 # 2 + -I 3 ^cos 0 + — figd cos 0 (5 9 2) 


5 - 20 



Hamilton’s Principle writes: 

t r l S 

f> J Qli<£ 2 sin 2 0 + ^hO 2 + - pgdcos^dt = [p 5 i] (593) 

t; <, 

Variations with respect to 8 , <j>, and rp yield to three equations of motions. 

88 : — I i4» 2 sin 28 — 1^8 — sin 0^ cost? -f ifj -+- f*gd cos 8 = 0 

8<f> : I\<j> sin 2 8 + Iz(j> cos 8 + ipj cos 8 = P# = Constant (5.9.4) 

6il> : J3 cos 8 + tpj = = Constant 

These equations are solved by Gear’s method in IMSL routine. The error tolerance is 
assigned as 10“ 10 . The solutions obtained by this method will be called as “analytical 
solution”. 

Three sets of initial conditions were investigated and compared with another numer- 
ical solution obtained by FJ. Mello [5-1]. The initial conditions are chosen to yield 
three different types of motions. 

Case 1 exhibits precession which is always in the same direction throughout in the 
motion; Case 2 exhibits precession which changes sign during the motion; and Case 3 
exhibits precession which does not reverse its direction but it does stop at points in its 
motion (cuspidal motion). 

The input data for each case is the same as Ref [5-1] given as: 

Case 1: 

The mass is 1.0. The axial moment of inertia is .40 

The transverse moment of inertia is .75. 

The initial orientation is in the yz plane +10 degrees from vertical. 

The initial angular velocity is (0., .9888, 7.5167) rad/sec. 

The distance from the mass center to the support is .2 

Gravity is 3.0 

Case 2: 

The initial angular velocity is (0., .20905, 6.2964) rad/sec. 

All other data are the same. 

Case 3: 

Initial angular velocity is (0., 0., 6.3794) rad/sec. 

All other data are the same. 

In this work the constraints are imposed via a Lagrange multiplier technique, and 
uses a fictitious stiffnesses (5.3.8, 5.4.6). The results are sensitive to the choice of 
these Lagrange multipliers. a\ is the fictitious stiffness associated with the normality 
condition of the Euler parameters representing the rigid body motion, and 02 that the 
hinge constraint. Several values were in uses as shown in table 5.1. 
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Table 5.1 Choices of fictitious stiffness 



Fictitious stiffness for 
Euler Parameters(ai) 

Fictitious stiffness for 
Hinge Constraint^) 

Case (A) 

10* 

10 7 

Case (B) 

10 6 

10 8 

Case (C) 

10 8 

10 10 


The motion of the spinning top was analyzed with 2,3, and 4 noded elements in 
time for the various motion types, and for the different values fictitious stiffness. The 
accuracy of the results is assessed by comparing the analytical solution with the numerical 
predictions for the x, y, and z coordinates of the center of mass of the spinning top. Table 
5.2 presents the results for the case 1 type motion. 100 time steps of At=0.06 sec were 
performed. For the two noded time element the D-hinge solution diverged for all values 
of the fictitious stiffnesses, the V-hinge solution was stable, but its accuracy was poor 
when compared to the results of ref [5-1]. For the three noded element, the D-hinge gave 
stable solutions but its accuracy is rather poor when compare to that of V-hinge which 
gave answers comparable to that of ref [5-1]. Finally the four noded element gave poor 
answers for the D-hinge, and very accurate predictions for the V-hinge. 

These results also clearly show the need to select appropriate values of the fictitious 
stiffnesses: high enough values should be selected as to properly enforce the constraint. 
Fig 5.1 shows the motion of the center of mass projected in the xy plane, for the two 
noded element It is that the numerical results are ahead of the analytical solution, in 
other words, the top “spins faster” in the numerical model. Fig 5.2 shows the result for 
the D-hinge, however, at Af=0.02 sec was selected, as Af=0.06 sec yield unstable results. 

Fig 5.3 and 5.4 present the results of the three and four noded elements with time 
increment. Fig 5.5 through 5.9 present the results for case 2 and 3 motions. It is clear that 
the D-hinge shows oscillations of increasing amplitude, typical of numerically unstable 
behavior, whereas the V-hinge yields a numerically stable solution. 

Fig 5.10 to 5.12 show the distance between the actual position of the base of the 
spinning top, and the position it is constrained to be. Surprisingly, the V-hinge results 
maintains this distance to a zero value, as it should, whereas the D-hinge results oscillate 
about the zero value and eventually become unstable. Fig 5.14 through 17 show similar 
results for the case 2 and 3 motion types. 
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Table 5.2 Relative Displacement errors of The Top Case 1 


2-noded Element 


Relative Errors (%) 

X 

Y 

Z 

Ref [1] 

0.976 

1.345 

0.0539 

V-Hinge 

(A) 

88.90 

6.636 

2.844 

(B) 

129.3 

-7.130 

2.692 

(C) 

130.0 

-9.196 

2.672 


3-noded Element 


Relative Errors (%) 

X 

Y 

Z 

Ref [1] 

0.308 

0.176 

0.0068 

D-Hinge 

(A) 

44.96 

-10.15 

2.947 

(B) 

10.91 

4.687 

1.611 

(C) 

6.133 

5.235 

1.345 

V-Hinge 

(A) 

38.06 

14.98 

0.056 

(B) 

4.422 

-1.221 

0.040 

(C) 

0.350 

-0.103 

-0.0048 


4-noded element 


Relative Errors (%) 

X 

Y 

Z 

Ref [1] 

0.068 

0.095 

0.0022 

D-Hinge 

(A) 

- 27.36 

21.75 

-12.46 

(B) 

-77.03 

35.66 

-15.45 

(C) 

-76.51 

34.22 

-14.94 

V-Hinge 

(A) 

38.05 

14.96 

1.0025 

(B) 

4.124 

-1.129 

0.045 

(C) 

0.048 

-0.017 

0.0005 
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Section 5.9.2 Pendulum Problem 

The characteristics of the rigid link constraint will be assessed by studying the single 
and double pendulum problems. Rigid link imposes the constraint of a fixed distance 
between two points. For instance the pendulum problem can be seen as the motion of 
a point mass in a two dimensional space subjected to the constraint of a given distance 
between the origin and the mass point. In the D-link the constraint condition is that of 
this fixed distance whereas in the V-link the constraint condition is the orthogonality of 
the position and velocity vector for the mass point. 

The analytical solution for the single and double pendulum are readily obtained and 
integrated. Three problems were analyzed with the following characteristics. 

Problem I : fig= 1.0 N, V| = 1.0 m/sec, L = 0.2m 

Problem II : pg = 10.0 N, V , = 1.0 m/sec, L = 0.2m 

Problem HI: fiig = ^2<7 = l-0 N, Vu = 1.0 m/sec,Vi 2 = 0.0 m/sec, 

L 1 = 0.2 m, 1 / 2 = 0 . 1 m. 

where ng is the gravity forces, V t the initial velocity, and L the length of the 
pendulum. 

Table 5.3 Choices of fictitious stiffness 



Fictitious stiffness for 
Velocity Link(a w ) 

Fictitious stiffness for 
Displacement Linkfc^) 

Case (A) 

10 5 

10 10 

Case (B) 

10 7 

10 12 

Case (C) 

10 10 

10 14 


Table 5.4 Relative Displacement errors 
Problem I 



X 

Y 

D-Link 

(A) 

58.79 

22.46 

(B) 

19.06 

38.67 

(C) 

18.53 

37.15 

V-Link 

(A) 

2.61 

4.60 

(B) 

0.40 

0.51 

(C) 

0.43 

-0.56 
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Problem II 


Relative Errors (%) 

X 

Y 

D-Link 

(A) 

*** 

*** 

(B) 

22.62 

15.13 

(C) 

21.95 

14.62 

V-Link 

(A) 

2.27 

4.95 

(B) 

0.28 

0.09 

(C) 

0.30 

0.15 


Problem III 


Relative Errors (%) 

X 

Y 

D-Link 

(A) 

*** 

*** 

(B) 

25.48 

5.50 

(C) 

23.18 

5.08 

V-Link 

(A) 

23.84 

0.95 

(B) 

1.27 

0.28 

(C) 

1.53 

0.30 
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Table 5.3 shows the three different cases of fictitious stiffnesses. From a number of 
trial and error reasonable ranges of fictitious stiffnesses are obtained. For D-link solutions 
for the single pendulum are going to be unstable in the case of a d < 10 9 . In the range of 
10 10 - 10 14 the solutions of D-link are converging for the single pendulum even though 
not converging for the double pendulum. But all the solutions are not accurate. On the 
other hand V-link does not seem to be sensitive to the values of fictitious stiffness a v if 
a v is large enough. The fictitious stiffness a v is investigated in the range of 10 5 - 10 10 . 

Table 5.4 shows the accuracies of the displacements for 2 different links. Again 
V-link has much more accurate solutions than D-link. In all 3 cases V-link has solutions 
with a reasonable accuracy. For problem I and II which is for the single pendulum 
D-link of case (C) yields the convergent solutions though not accurate. But for problem 
IE which is for the double pendulum D-link cannot obtain the convergent solution while 
V-link (B) and (C) has a quite accurate solution. 

Fig 5.18 — 5.20 show the results of problem I. Fig 5.17 shows the displacements 
of the pendulum along the time history with D-link. From fig 5.18 the solutions are 
quite sensitive to the values of a d . All three cases the solutions are oscillating. The best 
result is from case(C) the accuracy of which is at most 18%. Fig 5.19 shows the results 
with V-link. All three results are quite accurate and the best results are from case(B) 
the accuracy of which is as good as 0.5%. Those best results from D-link and V-link 
are investigated from now on. Fig 5.20 is the best result of problem I from D-link and 
V-link. Even though it looks identical the difference of the accuracies between the two 
links is large as seen from the table 5.4. Fig 5.21 shows the results of problem H. After 
200 time steps the pendulum oscillates for 5 periods. Fig 5.22 and 5.23 show the results 
of problem HI which is for double pendulum. From fig 5.22 displacements of D-link are 
diverging after 165 time steps in all three cases. With this A f =0.01 secs no convergent 
solution was obtained. Fig 5.23 shows the results with V-link which are quite accurate 
even after 400 time steps. In this double pendulum problem a„ = 10 s is not large enough 
to obtain the accurate solutions. But once a v > 10 7 the accuracies are within 1.5 %. 

Fig 5.24 shows how well the link preserves the distance from the mass point from the 
origin. Again the length of V-link remains the same as accurate as 10~ 5 error while the 
length of D-link is oscillating in the range of 5 % for problem I and 3 % for problem II. 
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Fig 5.1 Spinning Top Case 1: AT = .06 sec 
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Fig 5.2 Spinning Top Case 1: AT = .02 sec 
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Fig 5.3 Spinning Top Case 1: AT = .06 sec 
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Fig 5.4 Spinning Top Case 1: AT = .06 sec 
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Fig 5.5 Spinning Top Case 1: AT = .06 sec 
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Fig 5.6 Spinning Top Case 2 
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Fig 5.8 Spinning Top Case 3 
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Fig 5.10 Relative Errors at The Hinge(Case 1) 
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Fig 5.11 Relative Errors at The Hinge(Case 2) 
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Fig 5.12 Relative Errors at The Hinge(Case 3) 
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Fig 5.13 Spinning Top Case 2: AT = .06 sec 
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X Displacement * 10e-01 


Fig 5.14 Spinning Top Case 2: AT = .2 sec 


5-40 



Y Displacement * 10e-01 

- 0.5 - 0.4 - 0.3 - 0.2 - 0.1 0.0 0.1 0.2 0.3 0.4 0.5 



- 0.5 - 0.4 - 0.3 - 0.2 - 0.1 0.0 0.1 0.2 0.3 0.4 0.5 


X Displacement * 10e-01 

Fig 5.15 Spinning Top Case 2: AT = .2 sec 
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Fig 5.17 Spinning Top Case 3: AT = .1 sec 
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Fig 5.18 Disp. Link : Vi(x) = 1.0 , /xg = -1.0 
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Fig 5.19 Vel. Link : Vi(x) = 1.0 , yu.g = —1.0 
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Fig 5.20 Rigid Link : Vi(x) = 1.0 , fxg = -1.0 
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Fig 5.21 Rigid Link : Vi(x) = 1.0 , yug = -10. 
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Fig 5.22 Disp. Link : Vil(x) = 1.0 , fxg = -1.0 
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Fig 5.23 Vel. Link : Vil(x) = 1.0 , jxg = -1.0 
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Chapter 6 


Full Finite Element Modeling 


Section 6.1: Introduction 

The previous chapters have described the formulation of multibody dynamic problems 
with a finite element based modal analysis procedure to model each elastic body. This 
kind of formulation is very cumbersome, as demonstrated by the lengthy derivations and 
expressions of the last chapters. Furthermore, there is a tremendous overhead associated 
with modal methods consisting of the management and manipulation of all the coefficients 
of the modal expansion. This results in an increased computational cost which burdens 
modal analysis. 

An alternate formulation is to deal directly with a full finite element model, such as 
that used in chapter 4 for determining the reference solution. In this case, there is no 
need to distinguish between rigid body motion and elastic motion: the actual motion of 
each node is tracked by the finite procedure. Large rigid body motion results in finite 
rotations, however, finite rotations were already required to properly model the large 
elastic rotations. In this effort the Milenkovic parameters were used to model the finite 
rotations (see appendix A). This representation of the finite rotations is preferable to the 
Euler parameters as it involves three parameters only as opposed to four. This results 
in improved computational efficiency. 

An additional advantage of the full finite element formulation is that constraint 
equations (“hinges”) between two elastic bodies are much easier to formulate. In deed 
the constraint only involves the degrees of freedom of two nodes, on of each body. This 
contrasts with the formulation of constraint equations between elastic bodies modeled 
with a modal representation which involves all the elastic and rigid degrees of freedom 
of both bodies. Section 2 describes the formulation of a typical hinge in the full finite 
element model. 

Section 6.2: Hinge Element 

Consider a hinge with two components which undergoes a relative rotation. The two 
components of the hinge are defined by the triads e A and e B , respectively, which are 
coincident before deformation, i.e. e A = ej®. After deformation the triad e A has become 
e^ A (rotation tensor R(a)), and the triad e B has become e^ B (rotation tensor R(b)). The 
relative rotation <f>, between the two components is shown on Figure 6.1. The finite 
rotation tensor R(t) defines the initial orientation of the hinge: 

e ni = e n» = Rijifyinj = Cnj (6.2.1) 

and after deformation, we have: 

e ni Rij{ a ) = e n j j C n ,' = Rij(b) = C«j (6.2.2) 
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Fig 6.1 Relative Rotation at a Hinge 


Since the hinge only allows a relative rotation about e A = ef, the rotation tensor 
R(a) and R(b) are related by the following constraints: 


Ci = e[ B .e^ A = 0 
Ci = e 2 B .e^ A = 0 


(6.2.3) 


If the value of the relative rotation is desired, it can be obtained from the following 
constraints: 

C3 = 4 tan -1 ^ (?■ e* A ^j - <f> = 0 (6.2.4) 

where r is the relative rotation vector, and R(r) the corresponding finite rotation tensor. 
From (6.2.2), it is clear that 


<! = R ti (b)R kj (b)e-J = ^(rje^ 
or 

Rij(r) = R ik (b)R jk (a ) 

This tensor relationship can be expressed in the e* n A triad, as: 

R(rW) = R(b^R T (a^ 


(6.2.5) 

( 6 . 2 . 6 ) 

(6.2.7) 
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The right hand side can be transformed to the i n triad to yield 

J?(rW) = R T (t)R T (a)R(b)R(t) (6.2.8) 

where a,b,t, are the components of the rotation vector in its base triad i n . Using the 
formula for composition of rotation, (Appendix A), the Milenkovic parameter r of the 
relative rotations are found as 

r = jjR T (t) (aob - boa + baj;D = ao&o + a T b + (4 - ao)(4 - bo) (6.2.9) 
The constraints are now evaluated to yield: 

Ci = ej R T {b)R{a)e A = 0 

C *2 = e 2 R T (b)R(a)e^ = 0 (6.2.10) 

Co = 4 tan -1 — — bog? + a T b T ^e^ — <f> = 0 

It is preferable to enforce the two first constraints in a differential form, whereas the last 
constraint can be enforced as is since it is only used to define <j > . 


C h T 2< ^ T 
Ci = a 

4 - ao 
C 2 = a T 2 ° T 


:T 2 G T 


e;ei — b - — eoe, = 0 

4 - bo * 


■ejCj — £ 


:T 2 G T ~ 


= 0 


4 - ao 4 - &o 

Co = 4 tan -1 — — bog? + g?b T ^je^ — <f> = 0 

where e* = ^(6)^ , e£ = R(b)e 2 , e$ = i?(a)e^ 

These constraints could be enforced using a penalty technique: 

if 


( 6 . 2 . 11 ) 


— J 2 ( ai ^l + a2 ^2 + “3^3 ) + 


df 


( 6 . 2 . 12 ) 


where ai, 0-2, and ao are the penalty coefficients, and k is the stiffness of a torsional spring 
that can be present in the hinge. It is preferable, however, to use a mixed formulation 
to enforce the constraint (6.2.12) then transforms to 

if 


" / H 

u 


F\C\ + F 2 C 2 + F 3C3 + T<{> 




2\a, 


H h 

2 a.3 


2 it 


dt (6.2.13) 


where F\,F^ and F3 are the Lagrange multipliers associated with the three constraints, 
and T the torque in the spring. 
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Chapter 7 


Conclusions and recommendations 

Two approaches have been developed to analyze the dynamic behavior of multibody 
systems. In the first approach, each elastic body is represented in a local, noninertial 
frame of reference with unknown rigid body motions with respect to an inertial frame, 
and a finite element based modal analysis methodology is used to model the elastic 
behavior of the body. Constraint equations are used to model the interaction among 
the various elastic bodies. In the second approach, all elastic bodies are represented 
directly in a single inertial frame, and a full finite element methodology ( without a 
modal reduction) is applied. Constraint equations are used again to model the interaction 
among the elastic bodies. 

The following conclusions can be drawn from the study of the finite element based 
modal analysis methodology: 

1. The accuracy of modal methods strongly depends on the choice of the modal basis. 

2. Nonlinear kinematic couplings are poorly represented by natural vibration mode 
shapes. This is easily understood since both phenomena are of a different physical nature: 
one is a purely nonlinear kinematic phenomenon, the other a purely linear vibratory 
phenomenon. Even a large number of orthogonal vibration modes do not “synthesize” 
properly the nonlinear kinematic behavior. 

3. Adding perturbation modes to the classical natural mode shapes considerably im- 
proves the accuracy of modal methods. Perturbation modes contain information about the 
nonlinear behavior of structures extracted from higher order derivation of the Lagrangian. 

4. The nonlinearities associated with rotational dynamic effects are sometimes poorly 
respected by both natural vibration mode shapes and perturbation modes, resulting in a 
poor correlation for the angle of attack. 

5. When accurate predictions of rotor behavior is sought, modal analysis should be 
avoided, and full finite element methods should be preferred. 

6. A tremendous amount of overhead is associated with nonlinear modal methods. This 
involves the storage and manipulation of the many coefficients appearing in the elastic 
modes. The number of coefficients grows as N n where N is the number of modes, and 
n the highest power of the nonlinearities. 

7. When rigid body motions are added to the elastic behavior this overhead increases 
roughly tenfolds. This tremendous overhead is responsible for the very rapid increase 
in computational effort required to deal with modal methods as the number of modes 
increases. 

8. The computational effort involved in the integration of the full finite element calcu- 
lations presented in this work does not seem to be prohibitive when compared to that of 
the modal analysis. This observation should not be generalized. It is clear that as the 
number of degrees of freedom in the finite element model increases the cost of solution 
increases as well and will eventually become more expensive than that of the modal 
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solution. It seems however, that for typical rotorcraft problems the full finite element 
method is directly competitive with modal solutions. 

Two approaches were developed for dealing with the kinematic constraints. In the 
first approach the kinematic constraints are enforced as is, whereas in the second approach 
the time derivative of the constraints are enforced. In both cases a Lagrange multiplier 
technique is used to enforce the constraints within the framework of a mixed formulation. 

The following conclusions can be drawn from the study of the constraint equations 
enforcement: 

1. Enforcing the time derivative of the kinematic constraints yields numerical schemes 
which are far more accurate than those obtained from enforcing the kinematic constraint 
itself. 

2. When kinematic constraints are enforced, the problem becomes very “stiff’ due to 
the presence of the larger fictitious stiffness associated with the constraint. Integration 
schemes applied to these very stiff problems can easily become unstable, and this behavior 
was observed in the various examples treated here. When the time derivative of the 
constraint was enforced, this numerically unstable behavior disappeared. 

3. Enforcing the time derivative of the kinematic constraints is only slightly more 
complex than enforcing the kinematic constraint itself. 

The two methods developed in this study for the analysis of multibody dynamic 
systems differ only by their modeling approach for a single elastic body. The first 
approach relies on a modal approximation, the other on a full finite element model. 
The use of a modal approach requires modeling each elastic body in a local coordinate 
system. A large overhead is associated with nonlinear modal analysis, which puts this 
approach to a disadvantage when high order nonlinearities are present, and specially in 
the context of multibody analysis. Furthermore the accuracy of modal methods tend to 
deteriorate when nonlinearities are prominent. This discussion leads to the following 
recommendations for future work: 

1. For dynamic systems where nonlinear effects are not too pronounced, the finite 
element based modal analysis option can be pursued, and the following features of the 
modeling approach could be improved. The order of the nonlinearities could be reduced 
by ignoring higher order terms in the expression of the Lagrangian expression. For 
instance, keeping only quadratic terms would result in a linearized modeling of the system, 
keeping terms up to the third order would correspond to a “moderate rotation” type 
approximation. Such simplification would considerably decrease the overhead associated 
with the modal approach and make it much more computationally efficient at the expense 
of limiting its range of validity. 

2. Constraint equations are key to efficient multibody dynamic analysis. Enforcing the 
time derivative of theses constraints appears to improve numerical stability and accuracy. 
Further insight could be gained by computing the internal forces associated with the 
constraints. For instance, in the case of a rigid link, the load in the link is often 
a quantity of primary interest which can be obtained from the time derivative of the 
Lagrange multiplier used to enforce the constraint. 
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3. The full finite element approach should be used in complex multibody dynamic 
presenting strong nonlinearities. The resulting equations are rather “stiff”, even when 
the kinematic constraints are enforced in a time derivative fashion. To avoid numerical 
instabilities in the integration process, it is desirable, and probably necessary to use an 
integration scheme that provides high frequency numerical damping. Finite element in 
time procedures based on the time discontinuous Galerkin method might be well suited 
for this type of problems. 

4. The determination of the physical stability boundaries of a multibody system is an 
important problem. With a modal analysis involving a small number of degrees of 
freedom classical techniques such as the characteristic exponent method or the Roquet 
theory approach can be readily used. However, with the full finite element approach such 
method breaks down because of the presence of a large number of degrees of freedom. 
Innovative methods should be developed to deal with this situation. 
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Appendix A: Parametrization of Finite Rotations. 


Section A.1: The Geometric Representation of Finite Rotation 

Consider a finite rotation of magnitude <f> about an axis u (unit vector) and an arbitrary vector a. 
Let the rotation <f>, u bring this vector to b. From Figure A. 1 it is clear that: 

b = OC + CB =|| b || cos a u+ || b || sin a [s cos ^ + isin <f>] (Al.l) 

where the unit vectors s and t are defined as: 


_ _ (u x a) x u 

S = t x u = ^77 — — —rr, , 


u x a 



Figure A. 1 Change of Basis Viewed as a Single Rotation 

The fundamental property of rotation is to preserve length, hence: 

(u ■ a) =|| a || cosa =|| b || cos a; || it x a ||=|| a || sina =|| b || sin a (A1.3) 

Equation (Al.l) now becomes: 

b = (u • a)u + (u x a) x u cos <f> + (u x a) sin <f> (A1.4) 

For a unit vector u, the following relation can be shown to hold: 

(a • u)u = a + u x (u x a) (A1.5) 
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so that (A 1.4) becomes: 


b = a 4- sin <f>{u x a) + ( 1 — cos <f>)u x (u x a). 


(A1.6) 


Section A.2: The Rotation Tensor 

We wish to transform the above geometric notation into a tensor, index notation. For practical 
implementations it is necessary to work with tensor components in a particular system. The following 
notations will be used: 


Geometric notation 

Tensor notation 

Tensor component notation 

u 

Ui 

U 

U * V 

u i v l 

T 

U V 

U X V = —v X u 

Sij(u)vj = 

uv = —vu 


Sii is a skew symmetric tensor which components are: 


S ij ( U ) 


0 


U 2 

US 

0 

-Ul 

-U2 

Ul 

0 


This tensor has the following properties that are readily verified: 

Si, j( u ) = 

and Sik(u)Skj(u ) is a symmetric tensor. When u, is a unit vector, then 


Sik(u)Ski(u)Sij(u ) = - Sij(u ) 


(A2.1) 


(A2.2) 


(A2.3) 


In tensor notation (A 1.6) becomes: 

b, = R v a 3 (A2.4) 

where the rotation tensor R t) is given as: 

Rij(u) = Sij + sin^5, ; (u) + (1 — cos <f>)Sik(u)Skj(u) (A2.5) 

and 

Rij( u) = Rji(-u) (A2.5a) 
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An interesting expression can be found by expanding the trigonometric function in series and using 
(A2.3): 

R,j(r) = 8ij + Sij{r ) + ^S ik (r)S k j(r) + ^S tt (r)S w (r)Sy(r) + • • • hot (A2.6a) 

or 

Rij(r ) = exp[5,j(r)] (A2.6b) 

where the rotation vector r, is defined as r, = <f>u t . 

The following fundamental property of the rotation tensor can be readily verified: 

RikRjk = RkiRkj = Sij (A2.7) 


which implies: 


det(Rij) = 1 (A?.8) 

The eigenvalues A and corresponding eigenvector e, of the rotation tensor are: 


A = 1; 

A = cos <f) i \/ — 1 sin <f>\ 

e, = u, 

e, = a, db \/—lb t 

(A2.9) 

where a; is an arbitrary vector normal to u% 9 i.e. 



= 0 and 

bi = S tj (a)uj 

(A2.9a) 


Section A3: The Angular Velocity Vector 

A time derivative of (A2.7) yields: 

RikRjk = — ^ RikRjk ^ ■> 

which implies that this tensor is skew symmetric, hence we can write: 

Sij(w) = RikRjki 

where w; is the angular velocity vector. With the helps of (A2.5) we find: 

S tJ (u) = <f>Sij(u) + sin <f>Sij(u) + (1 - cos <f>)[Sik(ii)Sjk(u) - Sjk(u)Sik(u)], 
where the following identity was used: 

Sik(u)Ski(u)Sij(u ) = 0 . 


An alternate writing is: 


u>{ = <f)Ui + sin^it, + (1 — cos<f>)Si } (u)uj. 


(A3.1) 

(A3.2) 

(A3.3) 

(A3.4) 

(A3.5) 


A — 3 



Section A.4: The Virtual Rotation Vector. 


The virtual rotation vector can be defined by analogy to the angular velocity vector as: 


Sij(W '>) = bR tk R jk , 


(A4.1) 


where 8ifi t are the components of the virtual rotation vector. Note that there, is no vector rp t such that 
8(ip t ) is the virtual rotation. Taking a variation of A3.2 and a time derivative of A4.1 yields: 


Sij(8u) — SRikRjk + RikSRjk 

Sij (w) = 6RitR, t + SRitRjt , <A4 - 2> 

whicl\ after subtraction becomes: 

Sij(Soj) = Sij ($0) + RitSRjk ~ SRikKjk • (A4.3) 


In view of the orthogonality of the rotation tensor A2.7, this becomes: 

S.jCM = (Sifj + S t i(6xp)Sij(u) - Su(u>)Sij(8ip). (A4.4) 

The following identity: 


Si,{a)S h (b) - Si,(b)Sij(a ) = S tJ (S kl (a)b,) 


(A4.5) 


then yields: 

Sij(6u ) = Sij (W) + Sij(S u (Sif>)ui) (A4.6) 

and finally: 


8ui = 6ifii - S ik (u)8rp k . (A4.7) 

This important relationship relates the variation of the angular velocity vector to the virtual rotation 
vector and its derivatives. 

Section A.5: Change of Basis 

Consider a basis b l defined by the unit vectors and b% and an arbitrary vector a. Let R(u l ) 

be a rotation vector applies to each one of these vector to yield basis b 2 defined by and 

the vector a 2 , (A2.4) yields: 

a 2 = R,j(u l )a) (A5.1) 


A — 4 


(A5.2) 


This tensor relationship can be expressed in component form as: 

a 2 w = 

where the notation ()W is used to express the components of a tensor in the basis 6 1 . 




Figure A.2 Change of Basis 

It is clear that the components of a 1 in basis 6 1 are identical to the components of a 2 in the base 
6 2 , i.e.: 


hence (A4.2) becomes: 


a™ = a 2 ^ 


a 2 W = 


(A5.3) 


(A5.4) 


Since the starting vector a 1 is arbitrary this relationship holds for any vector, and more generally 
we can write the transformation law of the components of a vector v as: 

or (A5.5) 

It is also clear that the rotation vector has the same components in the two basis, hence: 

R(u l ^=R(u l ^j (A5.6) 


Consider now a second order tensor such as T i} = mb. j, where a, and bj are two arbitrary vectors. 
In component form we have: 


tM = d l Wi 


j[ 2] = a [%[2]T 


(A5.7) 


A — 5 



(A5.8) 


With the helps of (A4.5) the transformation law for second order tensor component is found as: 

T [2] = R t (u^TWRfu 1 ^ 

Finally, consider a set of basis 6 1 , b 2 ■ • • , 6” obtained by successive rotation R(u l ) , R[u 2 ) • • • , R(u n ) 
and let R(r ) be the rotation from base b 1 to base b n . Repetitive application of (A4.1) yields: 



fly(r) = Rit (u"- 1 ')Ru (u"" 2 ) • • • R, m (u 2 ) R m , (a 1 ) 

(A5.9) 

In component form, the following expressions hold: 



= R(r^ = i2^u n_1[B] )i?^u"" 2 W) ••• J R(u 2[n ^i2^ u l[n]^ j 

(A5.10) 

*( r[11 > 

| = J*(rM) = Jl(u»- 1 M)a(t«*- 2 M) -••i2( u 2 W)i?( u 1 [ 1 ]), 

(A5.ll) 

=R{ 

(rM) = (ti 1 W) 12 (w 2 ®) ...i2^u n - 2 [ n - 2 ^ J R^ n - 1 l n - 1 ]), 

(A5.12) 


where the last expression can be obtained from the first or second through repetitive use of the second 
order tensor component transformation law (A4.8). These various relationships can also be viewed as 
composition of rotations: the rotations u 1 , u 2 , ...,it" -2 , u n_1 are applied successively to yield a 
single composed rotation r. 

Section A.6: Euler Parameters 

The rotation tensor (A2.5) is expressed in terms of the rotation vector u and the magnitude of the 
rotation <f>. An alternate representation is in terms of the Euler Parameters defined as follows: 

<f> . <j> 

qo = cos qi = u, sin i = 1, 2, 3 (A6.1) 

to yield the rotation tensor (A2.5) as: 

RiM = 8ij + ZqoSijiq) + 2S ik (q)S k j(q) (A6.2) 

OT 

9o + 9i - ?2 “ 93 2(9192 - 9093) 2(9193 + 9092) 

R(q) = 2(9192 + 9093) 9o - 9? + 9? - 93 2(9293 - 9091) (A6.3) 

. 2(9193-9092) 2(9293 + 9091) 9o“9i“ 92 + 93_ 

All trigonometric functions have now disappeared from the rotation tensor. It is important to note 
the redundancy in the representation since four parameters are used when in fact only three parameters 
are required. This redundancy is clear when considering the definition A6.1 which yields the normality 
condition: 

9o + 9i + 9? + 9? = 1- (A6.4) 


A — 6 


The angular velocity tensor (A3.3) becomes: 


Si } (u) = 2qoSij(q) - 2 q 0 Sij(q) + 2[S ik (q)S jk (q) - S jk (q)S ik (u)] 
and the corresponding angular velocity vector is: 

w, = 2qoqi - 2qoq, + 2S tJ (q)q } 
or 


w(q) = 2 


-9190 + 9091 - 9392 + 9293 

-9290 + 9391 + 9092 - 9193 
-9390 - 9291 + 9192 + 9093 J 


= 2 H(q)q 


where 



-91 

90 

-93 

92 

H(q) = 

“92 

93 

90 

-91 


.-93 

-92 

91 

90 . 


(A6.5) 


(A6.6) 


(A6.7) 


(A6.8) 


The components of the angular velocity vector in the rotating system are found by using (A5.5): 


a*(9) = R t (iM<i) = 2 G(,)j 


where 



-91 

90 

93 

-92 

G(q) = 

-92 

-93 

90 

91 


.-93 

92 

-91 

90 . 


(A6.9) 


(A6.10) 


These matrices present the following remarkable properties: 

H(q)H T (q) = G(q)G T (q) = /, 

R(q) = H(q)G T (q ); G(q) = R T (q)H( q y, H(q) = R(q)G(q)-, 
H T (q)H(q) = G T (q)G(q) = I, - qcf; H(q)g = G(q)q = 0 
R(q) = 2H(q)G T (q) = 2 H{q)G T (q) 

* = R(i)R T ii) = 2 H(q)H T ( q y, u- = R T (q)R(q) = 2 G(,)G T (q); 


(A6.ll) 

(A6.12) 

(A6.13) 

(A6.14) 

(A6.15) 


Section A.7: Composition of Rotations with the Euler Parameters 


Since all trigonometric functions have been eliminated from the rotation tensor and angular velocity 
vector when expressed in terms of Euler Parameters, all finite rotation operations can be expressed in 


A — 7 



component form using a purely algebraic formalism. This task is eased when the following matrices 
are defined: 


A{q) = 


B(q) = 


C{q) = 


D(q) = 


90 

-91 

-92 

-93 

91 

90 

-93 

92 

92 

93 

90 

-91 

.93 

-92 

91 

90 . 

90 

-91 

-92 

- 93 ' 

91 

90 

93 

-92 

92 

-93 

90 

91 

.93 

92 

-91 

90 . 

90 

91 

92 

93 " 

91 

-90 

93 

-92 

92 

-93 

-90 

91 

.93 

92 

-91 

- 90 . 


(A7.1) 


(A7.2) 


1 0 
0 i*(<z)J 


(A7.3) 


(A7.4) 


Let q T = ( qo , gi, <72, 93) and r T = (ro, n, r2, n) be the Euler Parameters of two rotations in specific 
axes. The following formulae are readily verified: 

G(q)r = -G(r) £ ; tf( 9 )r = -JJ(r) £ ; 

A(q)B T (r ) = B T (r)A(9) 


A( g )r = 5(r) £ ; A r ( g )r = C T (r) £ ; B T (9)r = C(r) £ ; 


D(q ) = A( 9 )B T (g) = 2? r (g)A( £ ) = C(g)C( 9 ). 

Furthermore, the normality condition results in the orthogonality of these matrices 

A T (q)A(q) = A(q)A T (q) = If B T (q)B(q) = B(q)B T (q) = If 
C T (q)C(q) = C(q)C T (q) = If D T (,)D(q) = D(q)D T (q) = If 

From the definition of the angular velocity vector we find: 

u = 2B T {q)qj A(u>) = 2A(q)A T (q) ] B{u) = 2B T (q)B(q) 
uA=2A T (q)q; A(u*) = 2A T (q)A(q) ] B(u>*) = 2B(q)B T (q) 


(A7.5) 

(A7.6) 

(A7.7) 

(A7.8) 


(A7.9) 


With the help of the above relationships, the composition of rotation formula A5.12 can be shown 
to imply: 


A(rW) = a(,*w) a(, 2 I 2 1) ■ ■ .,<(,-11-11) 

B(rW) = B(,-il-'l) • ••£(,«) b(,iW) 


(A7.10) 


A — 8 



or 


rM = l- 2 l)j»-l["-l] 

rW = b( 9 "- 1 I"- 1 ))b( 9 "- 2 I"- 2 ]) • • ■ s( 9 2|2| j 9 1 i 1 i 

(A7.ll) 

An alternative representation of composition of rotation such as A5.ll implies: 


A^r^ = • • • A^g 2 ^ 

B (r [1] ) = B (V [1] ) B (? 2 [1] ) • • • B (g n ~ 2[1] ) B (? n " lll] ) 

(A7.12) 

or 

r(« = a( ? ”- 1 I 1 ])^( 9 ”- 2 W) 

rW = B(V [1| )b(V 11) ) ...b(»- j[1 ]) 4 ~>P> 

(A7.13) 

Section A.8: Rodrigues’ Parameters 


Rodrigues’ Parameters can be defined in relation to the Euler Parameters as: 


r, = 2—, 
90 

(A8.1) 

so that their geometric interpretation is: 


4> 

r, = 2 Ui tg-. 

(A8.2) 

It is clear that this representation presents an obvious singularity r, -* 00 when <f> — ► 
A8.1 can be inverted to yield: 

±7r. Relationship 

_ 1 _ r, 

90 (1 + r 2 /4) 1/2 ' ^ 2(1 + r 2 /4) 1/2 

(A8.3) 

where: 


r 2 = r J + r\ + r% 

(A8.4) 


The rotation tensor follows from A6.3and A8.3: 

j + '-f -r> ^+r 2 

a(r) = T+T 7 /! ^ + rs '-'i + 'i-'l T" r > (A8 - 5) 

'l f+n i-?-4 + ?J 

The components of the angular velocity vector are obtained from A6.7 and A8.3: 

t * i = Gr (A8.6) 


A — 9 



where 


G = 


1 

l + r 2 /4 


1 



~ T i T i 1 
1 

¥ 1 J 


(A8.7) 


Section A.9: Milenkovic Parameters 


The Milenkovic Parameters can be defined in relation to the Euler Parameters as: 

4?. 


a, = 


1 + ?o 


so that their geometric interpretation is: 


a, = 4 ui tan — . 

4 


i = 0, 1,2,3, 




(A9.1) 


(A9.2) 


Note that these parameters present no singularity within the range -n < <f> < n. Relationship A9.1 
can be inverted to yield: 


?« = 


i = 0,1, 2, 3. 


4 -oo 

The parameter ao can expressed in terms of the other three as: 

ao = (16 — a 2 )/8 

where 

a 2 = of + a\ + aj 

This representation involves 3 parameters only aj , 02, and 03, however ao is used in the 
formulae to simplify the notation. 

The rotation tensor follows from A6.3: 

a o 1 a i ~ a 2 ~ a 3 2(0102 — 0003) 2(0103 + 0002) 

2(0102+0003) a 2 — a 2 + 02 — a 2 2(0203 - aoai) 

. 2(0103 - aoa 2 ) 2(a 2 03 + a 0 ai) ajj - of - a% + ajj 


(A9.3) 

(A9.4) 

(A9.5) 

various 


R(a) = 


1 


(4 - ao)' 


The components of the angular velocity vector are obtained from A6.7 and A9.3: 

2 


LJ — 


4-o 0 


Gh 


where 


G = 


1 


ao + ^ T“°3 T + fl 2 

T + °3 Q 0 + ^ ^ - ai 


T-°2 T+°l «0 + ^ J 


4-ao 

This matrix has the following remarkable properties: 

GG = R-, GG t = G t G = I ; Ga = G T a = a 


(A9.6) 


(A9.7) 


(A9.8) 


(A9.9) 


A — 10 


The components of the angular velocity vector in the rotating system are then: 


u 


4 - ao 


G 1 a 


(A9.10) 


The matrix G is orthogonal, hence the following properties can be obtained: 

1 


GG t = 7 ; 7 = 


- 4 — ao 


-Ha 


G t G = 7 *; 7 * = 


1 


where 


H = 


~ 4-a 0 

1 -f f 

f 1 

"f V 1 




This matrix has the following properties: 

GH T = H- G T H = H T - Ha = a 


(A9.ll) 


(A9.12) 


(A9.13) 


Other matrix functions of the Milenkovic Parameters also play an important role. The following 
matrices are defined for an arbitrary vector b: 


D(b) = 


1 


D*(b) = 


4 - ao 
1 

4 - ao 
1 


H T V - —a.b T 


1 

4‘ 


(ab — a.b 

1 ’ 

v -- J 

4 — ao 


4 


-i(^> 


(A9.14a) 


4 - ao 
1 

4 — ao 


Hb - -a.b T 
4 


(ab — a. 6 ^) 

-, 1 f 

v — / 

1 

O 

O 




(A9.14b) 


and 


E(b) = 


1 


4 - a 0 
1 

4 - ao 


Hb T --a.b T ; 


— b — - (ab + a. b^^j 


4 - a 0 


-b — ^ (a.l? + i.a T — (a T .bj ij j 


(A9.15a) 


The off-diagonal terms of the D matrices are skew symmetric as can be seen from expanding their 
definition: 


D(b) = 


do —d% c?2 

d% do —d\ 
— ^2 d\ do 


do 

di 

1 

l- 9 * 

-1 

1 1 

-f 

¥ 

d 2 

1 

1 

O 

0 


-1 

- a -t 

ds. 


L-f 

* 

-1 



'61' 


62 


,& 3 _ 


(A9.16a) 


A — ll 



D'(i) = 


(A9.16b) 


' <*0 

-d$ 

d\ ' 


cIq 

d\ 

l 

\~ a i 

1 

i i 

i 

uo 

°o 

d\ 

ft 

J. 

i 

> 

d % 

A. 

4 — ao 

¥ 

L-f 

l 

a i: 

4 

1 


In contrast the off diagonal terms of the E matrices are not skew-symmetric 
These matrices enjoy remarkable properties: 


b i 

h 

b?, 


GD(b)G T = D{Gb)- G T D(b)G = D(G T by, 
RD(b)R T = D(Rb ); R T D(b)R = D^R T bj ; 


GE{b)G T = E{Gb)- 
RE(b)R T = E{Rb); 

GD*(b)G T = D*(Gb ); 
RD*(b)R T = D*(Rb); 
Furthermore, they are related as follows: 


G T E(b)G = E ( G T b)- 
R T E(b)R = E^R T bj; 

G T D*(b)G = D * (G T bj; 
R T D*(b)R = D*(R r b ); 


D{b)c = E T (c)b; D* ( b)c = E(c)b- 
GD(b) = E(b ) ; D{ b)G = E (G T bj ; 

G T D*(b ) = E T (b)- D\b)G T = E T (Gb)- 
E(b)c = D T (c)b 


where both b and c are arbitrary vectors. 

These matrices appear when taking derivative of expressions containing Milenkovic 
For instance, one can readily verify that: 



Finally, derivatives of these matrices become: 

S(D(b)c) = X(b, c)6a + E T (c)Sb + D{b)6c 


(A9.17a) 


(A9.17b) 

(A9.17c) 


(A9.18) 


Parameters. 


(A9.19) 


(A9.20) 


A — 12 



where 


X(b, c) = - \ lie + (b T c) I + ( D(b)c).a 2 

4 4 — ao t V / 


where 


S(E(b)c) = Y(b,c)6a + D*(c)Sb + E(b)6c 
Y(b, c) = ^ + (fc)l + (. E(b)c).a 2 


(A9.21) 


(A9.22) 


(A9.23) 


Composition of rotation is easily obtained from the corresponding relationships for the Euler 
Parameters. Let /, and g, be the Milenkovic Parameters of two successive rotations, and a, the 
parameters of the resulting rotation, then A7.4 yields: 


a _ MJ)g 
4 - ao (4 - /o)(4 - go) 

Computing ao from the first equation we obtain: 

4A(/)g 

~ (4 - /o)(4 - go) - fogo - hgi 


(A9.24) 


(A9.25) 


Other useful relationships can be readily obtained from the corresponding relationships on the Euler 
Parameters. 

Within the range of — ir < <f> < n the Milenkovic Parameters present no singularities since 

—4 < ai < +4 and 0 < ao < 2 , (A9.26) 

furthermore the orthogonality of the matrix G always yields a well defined angular velocity. However 
restricting the range of <f> would limit the range of admissible rotations. Hence, then definition of the 
Milenkovic Parameters is generalized to allow any magnitude of rotation: 


U , 7T ' 

a, =4tt, tg I - + k - 
4u t tgf^ + (k-l)^ 


k even 


k odd 


(A9.27) 


where — 7 r < <f> < n , and due to the following trigonometric identities: 

^( 4 + k l) =^( 4 )’ keven ; 
tg (4 + (* ” 1 k odd ; 


(A9.28) 


it is clear that the bounds A9.16 remain valid. If, at an instant during the motion of the structure, the 
rotation reaches, let say, the upper bound, i.e.: 


<t> = 7T + e €>0 


(A9.29) 


A — 13 



or, in terms of the parameters: 


a 2 > 16 (A9.30) 

one passes from the current state to a complementary one (i.e. from an even state to an odd one) to 
remain within bounds: 


<f> —> 4>' = <t> - 2 t; k—>k' = k + 1 (A9.31) 

The corrected value of the parameters is: 

a' = Aui tgf^ + (*' - 1)|) = 4 u, tgf^- ~ + k^\ 

_ A U n\ . (<t>\ 2 (A9 - 32 ^ 

— 4u, 2 J ~ ^ u */ — ~ 16aj/ a 

a' 0 can be shown to transform similarly so that: 

aj = -l 6a t /a 2 i = 0,1, 2, 3 (A9.33) 

and finally: 

o' = 16/a (A9.34) 

which clearly shows that the passage to a complementary state decreases the norm of the Milenkovic 
Parameters. 
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